{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4e4de91f-c4da-4bd1-852d-c01e83d69f0a",
   "metadata": {},
   "source": [
    "## Drone-Network Design Queueing Model\n",
    "\n",
    "This notebook uses the model and data from the following studies:\n",
    "\n",
    "- Lejeune M.A., Margot F. 2022. Drone-Network Design for Out-of-Hospital Cardiac Arrests. Working Paper, submitted.\n",
    "- Custodio.J., Lejeune M.A. 2022. Spatiotemporal Data Set of Out-of-Hospital Cardiac Arrests. *INFORMS Journal on Computing* 34 (1), 4-10.\n",
    "\n",
    "The data and its full description can be found here:\n",
    "https://github.com/INFORMSJoC/2020.1022\n",
    "\n",
    "The example reads from the VBOHCAR Excel file that contains data for:\n",
    "- out-of-hospital cardiac arrest (OHCA) incidents\n",
    "- potential base stations for drones\n",
    "- pairwise incident/base station (Harversine) distances\n",
    "\n",
    "The MILP reformulation is modeled in gurobipy and solved using the Gurobi Optimizer. \n",
    "\n",
    "To begin, we import the necessary packages for data intake, modeling, and visualization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a31dda-8e54-4d63-b109-c711b699f813",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "adfb071f-78d2-487c-8da4-99186c8c04b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import collections as mc\n",
    "import numpy as np\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "from datetime import datetime\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f97e2f4-0abd-4a54-a2bd-c2b94a9f44c0",
   "metadata": {},
   "source": [
    "### Input Data Processing and Visualization\n",
    "We use pandas to quickly read in the required sheets from Excel. There is a small amount of data prep needed, where some duplicate events need to be removed and text should be cleaned. For reference, an incident location encompasses a street block. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8e5103af-4fe3-4c99-a656-74ecf4e6eea4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read incident and base station data\n",
    "ohca_read = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/drone_network/OHCAs.csv', index_col=0)\n",
    "stations_read = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/drone_network/BaseStations.csv', index_col=0)\n",
    "\n",
    "# read data on indicent-base distance\n",
    "dist_read = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/drone_network/Base_to_OHCA_distance.csv', index_col=[0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bffa34fa-b159-483e-bc17-5cf31f68652b",
   "metadata": {},
   "source": [
    "The above reads the complete data set, including all OHCAs, base stations, and respective distances. We will subset indicent data to the last month and select seven stations to get an instance that will solve in the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "70a7760d-f53c-48e7-938a-52a6bade8866",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of days: 30 \n",
      "Number of OHCAs: 69\n"
     ]
    }
   ],
   "source": [
    "# select a subset of OHCAs so the model can be solved in this notebook (if you have a full license then you can use all incidents)\n",
    "start_date = '2019-05-31' \n",
    "end_date = '2019-06-30'\n",
    "n_days =  (datetime.strptime(end_date, \"%Y-%m-%d\") - datetime.strptime(start_date, \"%Y-%m-%d\")).days\n",
    "ohca_df = ohca_read[(ohca_read.ReceivedTime >= start_date) & (ohca_read.ReceivedTime <= end_date)]\n",
    "\n",
    "# select a subset of base stations so the model can be solved in this notebook (if you have a full license then you can use all stations)\n",
    "stations_df = stations_read.loc[[2,3,5,14,25,31,35]]\n",
    "# some bases are different types (e.g. fire and police), add the type to create a different base name\n",
    "stations_df['Base'] = stations_df.Street + ' ' + stations_df.Type\n",
    "\n",
    "dist = dist_read.loc[dist_read.index.get_level_values(1).isin(ohca_df.index)]\n",
    "dist = dist.loc[dist.index.get_level_values(0).isin(stations_df.index)]\n",
    "print(f\"Number of days: {n_days} \\nNumber of OHCAs: {ohca_df.shape[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9933c5dd-afd3-4025-b130-c2da5fb5ac85",
   "metadata": {},
   "source": [
    "Next, we create scatterplots to visualize the OHCA locations and stations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "02eb637f-34c8-4c9e-85ec-a4b15a02d480",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x26103255190>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnzUlEQVR4nO3da5BcZ33n8e9v7heNNaMrI0tCElgu2ypJhLGjpWKS2CYIVTYO2g0xKe86C7sOLi4bsiRAUbWEUFQZx1T2xe7CKqwrrkowaEE2FBeDILFZqoKUsS0JWbaxjW1p7FldLI2sGc19/vuiz0itnp6Znpnu6dvvU9U13c+59PN095z/eS7nOYoIzMzM0tUUOwNmZlZ6HBzMzGwKBwczM5vCwcHMzKZwcDAzsynqip2BfFixYkVs2LCh2NkwMysrTzzxxJmIWJltWUUEhw0bNtDd3V3sbJiZlRVJr0y3zM1KZmY2hYODmZlN4eBgZmZTODiYmdkUDg5mZjZFRYxWMrM0547AiX1w8Ti0rId1u6Fja7FzZWXGwcFsBr19gxzu6ePswAjLWhvYtradzvbmYmdreueOwDP3Q0MHNK+FkXOp19d9wgHC5sTNSmbT6O0bZP+xkwyOjLNiSSODI+PsP3aS3r7BYmdteif2pQJDQweo5vLzE/uKnTMrMw4OZtM43NNHW1MdbU311Ei0NdXT1lTH4Z6+YmdtehePQ/3SK9Pql6bSzebAzUpm0zg7MMKKJY2XXrdePMaG899DF1+BiRtKsy2/ZX2qKamh43La6PlUutkcuOZgNo1lrQ0MDI8BqcDw5tNfhuFzRFNaW/65I0XOZYZ1u1N5GzkHMXH5+brdxc6ZlRkHB7NpbFvbzoWhMS4MjbLy/Pe4GG0MRBurl7YUpC2/t2+QR4/28rUDr/Do0d759W10bE11Pjd0wGBP6q87o20e3KxkNo3O9mbedf1qDvf0oYuvEE1r2bS0hbam+tQKeWzLn+z8bmuqY8WSRgaGx9h/7CTvun713EdHdWx1MLAFc3Awm0Fne3Pq4DxxQ9KWX395YR7b8tM7v4FLfw/39JX20FmrWG5WMstFgdvyzw6M0Np45blaa2MdZwdG8rJ/s7lycDDLRYHb8tM7vycNDI+xrLUhL/s3mys3K5nlqoBt+dvWtrP/2EkgVWMYGB7jwtAYOzYtL8j7mc3GNQezEjDZ+d3cUMuZ/mGaG2rn1xltliez1hwkNQE/BRqT9b8ZEZ9Nln0U+AgwBnwvIv4iy/b/GfhPgIC/jYj/lqQvA74BbABeBt4XEeeSZZ8GPgiMAx+LiB8upJBm5eBS57dZCcilWWkYuCUi+iXVAz+T9AOgGbgd2BoRw5JWZW4oaQupwHATMAI8Kul7EfE88CngJxFxr6RPJa8/Kel64A7gBmAN8GNJmyNifOHFNTOzXMzarBQp/cnL+uQRwD3AvRExnKx3Ksvm1wE/j4iLETEGPA68N1l2O/Bg8vxB4PfT0r8eEcMR8RLwAqngYmZmiySnPgdJtZIOAaeA/RFxANgM3CzpgKTHJd2YZdOjwDslLZfUAuwC1iXLVkdEL0Dyd7LmcTVwIm0fPUlaZp7ultQtqfv06dO5FMPMzHKU02ilpElnu6R24OGkuagO6AB2ADcCeyVtiohI2+4ZSV8E9gP9wGFS/RMzUbYsZMnTHmAPQFdX15TlVvoOHz/H94/2cvKNIVZf1cSuLZ1sW98x+4ZmVnBzGq0UEX3AY8BOUmf0+5Jmp4PABLAiyzb/OyJ+LSLeCZwFnk8WnZTUCZD8nWyW6uFy7QJgLfDaXPJppe/w8XPs+elL9A+PsWZpM/3DY+z56UscPn6u2FkzM3IIDpJWJjUGJDUDtwHPAo8AtyTpm4EG4EyW7Vclf9cDu4GHkkXfAe5Knt8FfDst/Q5JjZI2AtcAB+deNCtl3z/aS3trHR0tjdTU1NDR0kh7ax3fP9pb7KyZGbk1K3UCD0qqJRVM9kbEdyU1AA9IOkpqJNJdERGS1gBfjYhdyfbfkrQcGAU+PDlcFbiXVFPUB4HjwB8ARMTTkvYCx0g1QX3YI5Uqz8k3hliz9Mphm0ub6nntfAnfZc2siswaHCLiCPC2LOkjwJ1Z0l8j1fE8+frmafb7OnDrNMu+AHxhtrxZ+Vp9VRPnh0bpaLl8M53zQ6OsvqqpiLkys0m+QtqKYteWTvoGxjh3cZiJiQnOXRymb2CMXVs6i501M8PBwYpk2/oO7n7nRpY01vHa+UGWNNZx9zs3erSSWYnwxHtWNNvWdzgYmJUo1xzMzGwKBwczM5vCwcHMzKZwcDAzsykcHMzMbAoHBzMzm8LBwczMpnBwMDOzKRwczMxsCgcHMzObwsHBzMym8NxKZmYF1Ns3yOGePs4OjLCstYFta9vpbG+efcMic83BzKxAevsG2X/sJIMj46xY0sjgyDj7j52kt6/0b2rl4GBmViCHe/poa6qjrameGom2pnramuo43NNX7KzNysHBzKxAzg6M0Np4Zet9a2MdZwdGipSj3Dk4mJkVyLLWBgaGx65IGxgeY1lrQ5FylDsHBzOzAtm2tp0LQ2NcGBplIoILQ6NcGBpj29r2YmdtVg4OZmYF0tnezLuuX01zQy1n+odpbqjlXdevLovRSrMOZZXUBPwUaEzW/2ZEfDZZ9lHgI8AY8L2I+Iss238c+I9AAL8A/kNEDEn6BnBtslo70BcR2yVtAJ4BnkuW/TwiPjTvEpqZFVFne3NZBINMuVznMAzcEhH9kuqBn0n6AdAM3A5sjYhhSasyN5R0NfAx4PqIGJS0F7gD+LuI+MO09b4EnE/b9MWI2D7vUpmZ2YLMGhwiIoD+5GV98gjgHuDeiBhO1js1w3s0SxoFWoDX0hdKEvA+4Jb5FMDMzPIvpz4HSbWSDgGngP0RcQDYDNws6YCkxyXdmLldRLwK3A8cB3qB8xHxo4zVbgZORsTzaWkbJT2V7PfmuRfLzAqtt2+QR4/28rUDr/Do0d6yuLDLcpdTcIiI8aSZZy1wk6QtpGoEHcAO4M+BvUkt4BJJHaSanjYCa4BWSXdm7P79wENpr3uB9RHxNuDPgK9JuiozT5LultQtqfv06dO5FMPM8qScr/y13MxptFJE9AGPATuBHmBfpBwEJoAVGZvcBrwUEacjYhTYB7xjcqGkOmA38I209xiOiNeT508AL5KqpWTmZU9EdEVE18qVK+dSDDNboHK+8tdyM2twkLRSUnvyvJnUAf9Z4BGSfgJJm4EG4EzG5seBHZJaklrFraRGIk26DXg2Inoy3q82eb4JuAb41XwKZ2aFUc5X/lpuchmt1Ak8mBywa4C9EfFdSQ3AA5KOAiPAXRERktYAX42IXRFxQNI3gSdJDXd9CtiTtu87uLJJCeCdwF9JGgPGgQ9FxNmFFNLM8mvyyt+2pvpLaeVy5a/lRqnBSOWtq6sruru7i50Nm0W5Tl1sU032ObQ11dHaWMfA8BgXhsbK5gIvS5H0RER0ZVvmK6RtUbgDs7KU85W/lhvf7McWRXoHJnDp7+GePh9QylSxr/x1TbSwHBxsUZwdGGHFksYr0lob6zjTP1ykHFkpyvWAn96stWJJIwPDY+w/dtK1lzxycLBF4Q5Mm022A/43u0+woq2RgCuChWuihec+B1sU5Tx1sS2OzGsnRscnePnsRX558sKUfioPpS081xxKXKW0q052YB7u6eNM/zDLWhvYsWl5WZbFCiOz6fHlMxdZ2lTP6MTEpQvtIBVEXBMtPAeHElZp7arF7sC00pZ5wL8wPEp9rWhrunyYmuyn+u1rV7H/2MlLaZNDaXdsWl6UvFciNyuVME9RYNUks+mxrkb0XRxjw/Ill9aZrB14KG3hueZQwjzCx6pJZtPj5tVLODOQqj1MREypHbgmWlgODiXM7apWbTIP+JN9bu6nWnwODiVs29r2vLarVkrntlUP1w6Kx30OJWymdtW53mjF01eY2Vy45lDisp05zWcUky8aMrO5cHAoQ/M50Gd2bp8dGOGlMxfoPT8E4CYmM7uCm5XK0HyuDp3s3J7c/tCJPt4YGuNNVzW7icnMpnBwKEPpB/pJs41iSh9D/tKZC0hBhNi0stXXT5jZFA4OZWg+8xSld273nh+irbGe7euWsqw11dTkeWnMLJ37HMrQfOcpSu/cHhwZ9/UTZjYtB4cytZDx3/m+fsLMKo+DQxUq5gypvhDPrDw4OFSpYlx5WmmzzJpVMgcHy6oQZ/i+EM+sfMw6WklSk6SDkg5LelrS59KWfVTSc0n6fdNs//Fk+VFJD0lqStL/UtKrkg4lj11p23xa0gvJvt+dj4Ja7go11Ybv3mXVYq7T25SiXGoOw8AtEdEvqR74maQfAM3A7cDWiBiWtCpzQ0lXAx8Dro+IQUl7gTuAv0tW+ZuIuD9jm+uTdW4A1gA/lrQ5IsbnV0Sbq0Kd4XuWWasG6c2nNRIHfvU63/9FL/9q0zJ+69ryaUKdNThERAD9ycv65BHAPcC9ETGcrHdqhvdoljQKtACvzfKWtwNfT/b7kqQXgJuAf54trzZ32ZqPCnUfCY+SKm/VMphgoeWcPLkaHQ+O9JynuaGGFa2N/PJkP6PjlE0fW04XwUmqlXQIOAXsj4gDwGbgZkkHJD0u6cbM7SLiVeB+4DjQC5yPiB+lrfIRSUckPSCpI0m7GjiRtk5PkpaZp7sldUvqPn36dC7FsAzTNR+JmPMV2Lnw3bvKV7XM6puPck42n778ej/NDTU019fR3FDL2ESU1UwEOXVIJ0062yW1Aw9L2pJs2wHsAG4E9kralNQ0AEgO+LcDG4E+4P9IujMi/h74MvB5UrWQzwNfAj4AKFsWsuRpD7AHoKura8pym910zUdDo+NcGEoFh3yf4Xt+/vJULYMJ8lHOyebTC0NjXHXpf2qCtsb6K2rgpV4Tm9P0GRHRBzwG7CR1Rr8vUg4CE8CKjE1uA16KiNMRMQrsA96R7OtkRIxHxATwt6Sajkj2uy5tH2uZvSnK5mG6DuIAn+HbFaplMEE+yjk5vU1djRgcHWNwZJzB0XE2rGi5VAMvh5pYLqOVViY1BiQ1kzrgPws8AtySpG8GGoAzGZsfB3ZIapEk4FbgmWSbzrT13gscTZ5/B7hDUqOkjcA1wMH5FM5mNtMEfp3tzezc0skf/fqb2bml04Ghys1nssdylI9yTjafbl7dxukLI0zEBFvXXkV9bc2lOdDSayg1UklOfplLs1In8KCkWlLBZG9EfFdSA/CApKPACHBXRISkNcBXI2JXRByQ9E3gSWAMeIqkKQi4T9J2Uk1GLwN/AhARTyejmo4l23zYI5UKYzE6iItSdT53BE7sg4vHoWU9rNsNHVsL+54VrloGE+SrnJ3tzbz/19/Mb1276tLvv7mh9tJMBP/03Ck21LzI6t7v0zT6KkP1V3Pyql28PPCWQhRrXpTWRVC2urq6oru7u9jZKEuFPHinD+lL/0craBPVuSPwzP3Q0AH1S2H0PIycg+s+4QCxQKXeRp4vi1HOn/3LT1h/8svQ0MFYTRt1Exdg5BzHV9/Db9x4a17fayaSnoiIrmzLfIV0lStEB/HkP9dPf3maxroaruu86lLVGQrciXliXyowNCSD3yb/ntjn4LBApTqYIN8H88Uo57b4R3qiDaKNBtUwEG0QY2yLfyTV+l58vp+D5VV6R1uNQIhDJ85zdiA1QqPgnZgXj6dqDOnql6bSreKUQ8duNm0Tvaxd9Sbqa8XFkTHqa8XaVW+ibaK32Fm7xDUHy2q+Z2PpHW1XNdczPDZBc30tL5+5yLLWxpw79+Z9NtiyPtWMNFljgFTTUsv62be1slO2Q2xb1tM2co62VWm/05Fz0FA6v1PXHGyKhZyNpQ8F3LB8CYMjEwQTvDE0mtMd6xb6/qzbnfonGzkHMXH5+brduRTdykzZDrHNw++00PM3OTjYFAsZZpc+FHBZawPb17UTAUHkfL3Egob5dWxNdT43dMBgT+pvgTujK2GStXJVtkNsF/g7XYzmNDcr2RQLmVspcyhgfa3YuGLJnEYoLXhup46ti9b57HtUFFdZD7FdwO90MZrTXHOwKRZyNpaP+ZPK6WywHC5mqmTVOl/XYjSnueZQJKU8ZnyhZ2MLHQpYTmeDhZrB1nJXqkNsC2kxpr93zaEISn34XbHPxor9/nNRTrUcqxyT8zddGBplIiLnwR5z4ZpDEZTD8Ltin40V+/1zVU61HKsckydQh3v6ONM/zLLWhktTc+SLg0MRuCmicizGP6lZNoU+gXJwKALfLrOylEstx2wu3OdQBIvRXmhmthAODkVQTh2uZlad3KxUJG6KMLNS5pqDmZlN4ZqDFUwpX+hnZjNzzcEKotQv9DOzmTk4WEF4ziGz8ubgYAVRtvPsmxng4GAF4jmHzMrbrMFBUpOkg5IOS3pa0ufSln1U0nNJ+n3TbP/xZPlRSQ9JakrS/1rSs5KOSHpYUnuSvkHSoKRDyeMreSqrLSJf6GdW3nKpOQwDt0TENmA7sFPSDkm/DdwObI2IG4D7MzeUdDXwMaArIrYAtcAdyeL9wJaI2Ar8Evh02qYvRsT25PGheZbNisgX+pmVt1mHskZEAP3Jy/rkEcA9wL0RMZysd2qG92iWNAq0AK8l6/8obZ2fA/92PgWw0uUL/czKV059DpJqJR0CTgH7I+IAsBm4WdIBSY9LujFzu4h4lVSN4jjQC5zPCAqTPgD8IO31RklPJfu9eZo83S2pW1L36dOncymG5ZHvm2xW2XIKDhExHhHbgbXATZK2kKoRdAA7gD8H9kpS+naSOkg1PW0E1gCtku7MWOczwBjwD0lSL7A+It4G/BnwNUlXZcnTnojoioiulStX5lpeywNfw2BW+eZ0hXRE9El6DNgJ9AD7kmang5ImgBVA+mn8bcBLEXEaQNI+4B3A3yev7wJ+F7g12Q9JM9VkU9UTkl4kVUvpnm8hLb/K4WZFZvlUjVf7zxocJK0ERpPA0EzqgP9FUv0QtwCPSdoMNABnMjY/DuyQ1AIMAreSHOQl7QQ+CfxmRFzMeL+zETEuaRNwDfCrhRXT8mnWmxWdOwIn9sHF49CyHtbtho6ti5a/avxHtsKZrCm3NdWxYkkjA8Nj7D92suIHWOTSrNQJ/JOkI8C/kOpz+C7wALBJ0lHg68BdERGS1kj6PkDSN/FN4EngF8n77Un2+9+BNmB/xpDVdwJHJB1Otv1QRJzNR2EtP2a8huHcEXjmfhg5B81rU3+fuT+Vvgjc5GX5Vq1X++cyWukI8LYs6SPAnVnSXwN2pb3+LPDZLOu9dZr3+xbwrdnyZYU109n3jPdNPv6/oKEj9YDLf0/sW5Tag5u8LN8KeVvfUq7l+gppm2K2s+8Zr2G4eBzql165w/qlqfRF4Gk7Cq/aRqoV6mr/Uq/lespumyKXs+9pr2FoWZ9qSpqsMQCMnk+lL0CuZ1gLuT93KZ/FLZpZ+ouqsf19xpryApR6Lbeqaw7VdgaUqwWdfa/bnQoOI+cgJi4/X7d73vmZyxnWfKftKNRZXFn9xnLoL6rG9vdCXe1f6rXcqg0OpV6lK6YFVaM7tsJ1n0jVHAZ7Un+v+8SC+hvmckCa7z9yIQ56ZfcbO7Hvcn+Rai4/P7Hv0iqlfkArlM72ZnZu6eSPfv3N7NzSmZcz+1KfnLJqm5VKvUpXTAuuRndszWvn81w7BOczbUchOh3L7jd28XiqxpAuo79oIc12dqVCNVflS9XWHKr1DCgXpTZp3mKcYRXiPcruN9ayPtU/lC6jv8iz7eZPqf2fZaramoPPgGZWSpPmLcYZViHeo+x+Y+t2p/oYIFVjGD2f6nd4ywcvrTJ5QDvc08eZ/mGWtTawY9PykvmtlJtS+j/LpGTWirLW1dUV3d1zm10jfdRF+sGglCK3XbYYI4ny/R5l+Rsr8tXttrgkPRERXVmXVWtwAA9dtMLzb8xK2UzBoWqbleBylW7yH/ifnjvlf+AyU+oH31JuNjCbSdV2SE8qu+GGdom/O7PCqfrgUI0X9VQKf3dmhVP1waHshhvaJf7uzAqn6oNDqV+laNPzd2dWOFUfHHxRT/nyd2dWOFUfHEr9KkWbnr87s8Kp6qGsk0pxuGGhhmiW+tDPuSrF786sElR9zaEUFXL6aA/9NLNcODiUoEIN0fTQTzPLlYNDCSrUEE0P/TSzXDk4lKBCDdH00E8zy9WswUFSk6SDkg5LelrS59KWfVTSc0n6fdNs//Fk+VFJD0lqStKXSdov6fnkb0faNp+W9EKy73fno6DlpFBDND3008xylUvNYRi4JSK2AduBnZJ2SPpt4HZga0TcANyfuaGkq4GPAV0RsQWoBe5IFn8K+ElEXAP8JHmNpOuTdW4AdgL/U1Lt/ItYfgo1RNNDP80sV7MOZY3UnN79ycv65BHAPcC9ETGcrHdqhvdoljQKtACvJem3A7+VPH8QeAz4ZJL+9WS/L0l6AbgJ+Oe5FKzcFWqIpod+mlkucupzkFQr6RBwCtgfEQeAzcDNkg5IelzSjZnbRcSrpGoUx4Fe4HxE/ChZvDoiepP1eoFVSfrVwIm03fQkaZl5ultSt6Tu06dP51KMktfbN8ijR3v52oFXePRor4eYmlnR5BQcImI8IrYDa4GbJG0hVSPoAHYAfw7slaT07ZJ+hNuBjcAaoFXSnbO8nbKkTbkjUUTsiYiuiOhauXJlLsUoab4Gobr4RMBK3ZxGK0VEH6nmn52kzuj3RcpBYAJYkbHJbcBLEXE6IkaBfcA7kmUnJXUCJH8nm6V6gHVp+1jL5aaoiuVrEKqHTwSsHOQyWmmlpPbkeTOpA/6zwCPALUn6ZqABOJOx+XFgh6SWpFZxK/BMsuw7wF3J87uAb6el3yGpUdJG4Brg4HwKV058DUL18ImAlYNc5lbqBB5MRgzVAHsj4ruSGoAHJB0FRoC7IiIkrQG+GhG7IuKApG8CTwJjwFPAnmS/95JqivogqSDyBwAR8bSkvcCxZJsPR8R43kpcoiavQWhrqr+U5msQKtPZgRFWLGm8Iq21sY4z/cNFypHZVEoNRipvXV1d0d3dXexsLMhkU0NbUx2tjXUMDI9xYWjMQ00r0KNHexkcGb/iRODC0CjNDbXs3NJZxJxZtZH0RER0ZVvmK6RLhK9BqHyTndC/Ot3PgZde58TZi74Y0UqWp+wuIYW8BqHQU3VX2lTg+TL5ufzqdD/Hz17k2tVX8dZVbTTV1/LcyQsMjY6zcWUrOzYt9+dlJcXBoQqkN1mtWNLIwPAY+4+dzFvNpND7L5aFBrz0z+WNwTFqa8Tzp/ppbaxj/bJWOloa3JRkJcvNSlVgPqNj5jIOvxJH3+RjuGn65zIwMkZ7cwPNDTW8/HpqwgGPRrNS5ppDFZjr6JjZagKZZ9QvnR7gLauW5Lz/cpB+YAcu/T3c05dz7SH9c29rrGdodIKm+lreGBoFPBrNZlbsplrXHKrAXKfqnqkmkO2M+pWzF3n13GDO+y8H+bjuJP1z37CihcHRcfoGR1jSWOtOaJtRKVwo6eBQBeY6VfdMB8ZsgWPz6iU8d/KNipoKPB/3vkj/3NtbGnjrqlbGJ+CqpgaPRrMZlUJTrZuVqsDkMNnDPX2c6R9mWWvDjKNjZrogL1sT1dqOFoZGxy8Nw51t/+Vg29p29h87CXDFdSc7Ni3PeR+Zn/ua9mbes6WzrD8XWxylcKGkg0OVmMsw2ZkOjId7+rIGjk0rl+R91E0x21znGlBn2o+Dgc1VKcyY4OBQoRZyYJ3twLjQM+pc8z+f4bH5DCg+sFux5KPmulCePqMCFXoqjsU4o5/PFBOegsQqyWL8n800fYZrDhUoH8MwZ7IYZ9TzaXMtdLnNFlOxa64ODhWoFDqzFmo+ba6VUG6zTMXqe/NQ1gqUj2GYxTbX4bdQGeU2S1fM6x0cHCrQfA6spWY+s9RWQrnN0hXzegc3K1WgfA3DLLa5trlWSrnNJhWzqdTBoUIVuzOrWKq13FaZinm9g5uVzMxKVDGbSh0czMxKVDHvEOlmJTOzElasplLXHMzMbIpZg4OkJkkHJR2W9LSkz6Ut+6ik55L0+7Jse62kQ2mPNyT9abLsG2npL0s6lKRvkDSYtuwr+SuumZnlIpdmpWHglojol1QP/EzSD4Bm4HZga0QMS1qVuWFEPAdsB5BUC7wKPJws+8PJ9SR9CTiftumLEbF9XiUyM7MFmzU4RGpmvv7kZX3yCOAe4N6IGE7WOzXLrm4lddB/JT1RkoD3AbfMLetmZlYoOfU5SKpNmn1OAfsj4gCwGbhZ0gFJj0u6cZbd3AE8lCX9ZuBkRDyflrZR0lPJfm+eJk93S+qW1H369OlcimFWUnr7Bnn0aC9fO/AKjx7tXdRbQJrNJqfRShExDmyX1A48LGlLsm0HsAO4EdgraVNkmQNcUgPwe8Cns+z+/VwZNHqB9RHxuqS3A49IuiEi3sjI0x5gD6Sm7M6lHGalYr73q7DSUMwbUS2WOQ1ljYg+SY8BO4EeYF8SDA5KmgBWANlO498DPBkRJ9MTJdUBu4G3p73HMKl+DiLiCUkvkqql+IYNaarhx1nJPL14+aqWwJ7LaKWVSY0BSc3AbcCzwCMk/QSSNgMNwJlpdpNZO5h0G/BsRPRkvF9t8nwTcA3wq9yKUx2KOVOj5cfZgRFaG688N2ttrOPswEiRcmS5KuZkeIspl5pDJ/BgcsCuAfZGxHeTpqIHJB0FRoC7IiIkrQG+GhG7ACS1AO8C/iTLvrP1Q7wT+CtJY8A48KGIODufwlUqn3WWv1K4R7DNT7XcNySX0UpHgLdlSR8B7syS/hqwK+31RSDrjU8j4o+zpH0L+NZs+apm1fLjrGSlcI9gm7vevkFeeX2Ap46fY2VbIxuWL7kU6CstsPsK6TLkm9qUv2LOmWPzM9mcu/qqJupqajg/OMpTx89x4uzFirxviOdWKkM+66wMnl48vwo1SGNyv//3+dM01NZwXedS3ra+nZfPXOR0/zD/741B/vgdGyvuu3TNoQz5rNPsSoUapJG+XyEkOHSiDxC/9uYO3nX9at68vLUi//dccyhTPus0u6xQgzTS93tVUz3DYxM0NwQvv97PstZlFd2c65qDmZW9Qg0NTt/vhhUtDI6OExPwxuBoxd+j3MHBzMpeoQZppO93WWsj29ctJQgmgopvznWzkpmVvUIN0sjcb31tDRtXLKnooDDJwcGsylXCVCyTgzQO9/Rxpn+YZa0N7Ni0fMHlKNR+y4GDg1kVq6R5ggo1SKNaB3+4z8GsilXLPEE2dw4OZlXMEwDadBwczKqYp2Kx6Tg4mFWxbWvbuTA0xoWhUSYiKn7svuXOwcGsinkqFpuORyuZVblqHY1jM3PNwczMpnBwMDOzKdysZFbtzh2BE/vg4nFoWQ/rdkPH1mLnyorMwcGsmp07As/cDw0d0LwWRs6lXl/3iYoKEJUwRchic7OSWTU7sS8VGBo6QDWXn5/YV+yc5U2hbgRU6RwczKrZxeNQv/TKtPqlqfQK4SlC5mfW4CCpSdJBSYclPS3pc2nLPirpuST9vizbXivpUNrjDUl/miz7S0mvpi3blbbdpyW9kOz73Xkqq5llalkPo+evTBs9n0qvEJ4iZH5y6XMYBm6JiH5J9cDPJP0AaAZuB7ZGxLCkVZkbRsRzwHYASbXAq8DDaav8TUTcn76NpOuBO4AbgDXAjyVtjojxOZfOzGa2bneqjwFSNYbR86l+h7d8sLj5yqPJKUImbx0KniIkF7PWHCKlP3lZnzwCuAe4NyKGk/VOzbKrW4EXI+KVWda7Hfh6RAxHxEvAC8BNs+XTzOahY2uq87mhAwZ7Un8rrDPaU4TMT06jlZKz/ieAtwL/IyIOSNoM3CzpC8AQ8ImI+JcZdnMH8FBG2kck/XugG/gvEXEOuBr4edo6PUlaZp7uBu4GWL++cqrAZouuY2tFBYNM1XzDnoXIKTgkTTrbJbUDD0vakmzbAewAbgT2StoUEZG5vaQG4PeAT6clfxn4PKlayOeBLwEfAJQtC1nytAfYA9DV1TVluZnNrJqGd5bKFCHl9JnPabRSRPQBjwE7SZ3R70uanQ4CE8CKaTZ9D/BkRJxM29fJiBiPiAngb7ncdNQDrEvbdi3w2lzyaWYz8/DOxVdun3kuo5VWJjUGJDUDtwHPAo8AtyTpm4EG4Mw0u3k/GU1KkjrTXr4XOJo8/w5wh6RGSRuBa4CDuRXHzHLh4Z2Lr9w+81yalTqBB5N+hxpgb0R8N2kqekDSUWAEuCsiQtIa4KsRsQtAUgvwLuBPMvZ7n6TtpJqMXp5cHhFPS9oLHAPGgA97pJJZfp0dGGHFksYr0lob6zjTP1ykHFW+cvvMZw0OEXEEeFuW9BHgzizprwG70l5fBJZnWe/fzfCeXwC+MFvezGx+PLxz8ZXbZ+4rpM2qkId3Lr5y+8wdHMyqkO8At/jK7TP3rKxmVapUhndWk3L6zF1zMDOzKRwczMxsCgcHMzObwsHBzMymcHAwM7MplGWevLIj6TQw21Tgi2UF008jUsmqtdzgsrvs5evNEbEy24KKCA6lRFJ3RHQVOx+LrVrLDS67y16Z3KxkZmZTODiYmdkUDg75t6fYGSiSai03uOzVqqLL7j4HMzObwjUHMzObwsHBzMymcHDIgaRvSDqUPF6WdCht2VZJ/yzpaUm/kNSUZfu/lvSspCOSHp687Wqy7NOSXpD0nKR3L06JcpeHsv9BsnxCUlda+gZJg2n7/soiFSlnhSp7sqxkv/c8lHuZpP2Snk/+diTp1fCdZy17sqxkv/OsIsKPOTyALwH/NXleBxwBtiWvlwO1Wbb5HaAuef5F4IvJ8+uBw0AjsBF4Mdv2pfKYZ9mvA64FHgO60tI3AEeLXaYilb1svvd5lvs+4FPJ80+l/d6r4Tufruxl851PPlxzmANJAt4HPJQk/Q5wJCIOA0TE65HlftcR8aOIGEte/hxYmzy/Hfh6RAxHxEvAC8BNhSzDfC2g7M9ExHOLl9P8K0DZy+J7n2+5SZXvweT5g8DvFzireVeAspfFd57OwWFubgZORsTzyevNQEj6oaQnJf1FDvv4APCD5PnVwIm0ZT1JWinKR9kzbZT0lKTHJd2cv6zmXb7LXi7f+3zLvToiegGSv6vSllX6dz5d2cvlO7/Ed4JLSPox8KYsiz4TEd9Onr+fy2cSkPr8fgO4EbgI/ETSExHxk2ne4zPAGPAPk0lZVlv0scWLUfYseoH1EfG6pLcDj0i6ISLemF8p5qdIZS/69+7vvPq+87lycEhExG0zLZdUB+wG3p6W3AM8HhFnknW+D/waMOUHI+ku4HeBWyNphEy2X5e22lrgtfmWYb4KXfZp3nMYGE6ePyHpRVJnZ91zLsACFKPslMD3XuByn5TUGRG9kjqBU8l7VsN3nrXslMB3PlduVsrdbcCzEdGTlvZDYKukluQH9ZvAscwNJe0EPgn8XkRcTFv0HeAOSY2SNgLXAAcLVoL5m3fZpyNppaTa5PkmUmX/VR7znC95Lzvl8b0vpNzfAe5Knt8FfBuq5jvPWnbK4zu/UrF7xMvlAfwd8KEs6XcCTwNHgfvS0r9KMkKFVOfTCeBQ8vhK2nqfITVy4TngPcUuZwHK/l5SZ03DwEngh0n6v0m2PQw8CfzrYpdzscpeDt/7Asu9nNQZ9fPJ32VV9J1nLXs5fOeZD0+fYWZmU7hZyczMpnBwMDOzKRwczMxsCgcHMzObwsHBzMymcHAwM7MpHBzMzGyK/w/LbYeLJbkdHQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(data = ohca_df, x = 'Longitude', y = 'Latitude', alpha = .3)\n",
    "plt.scatter(data = stations_df, x = 'Longitude', y = 'Latitude', alpha = 0.5, c = 'orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "60e12679-3154-47d6-bd8d-21e10e4d49c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a base / OHCA data frame using the location names\n",
    "df = dist.join(ohca_df[['Incident_Location']])\n",
    "df = df.join(stations_df['Base'])\n",
    "df = df.drop_duplicates()\n",
    "df = df.set_index(['Incident_Location','Base'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98dc111f-9a09-4deb-92fc-957e7c04fefe",
   "metadata": {},
   "source": [
    "### MILP Model\n",
    "\n",
    "#### Parameter Calculation\n",
    "We assume that drones at each station will have the same speed and mainenace times. The total time it takes a drone to service an OHCA is travel time (to and back from the OHCA location) plus maintenance time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d7e52800-a8f2-438c-8fb9-50adf5aee8db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# limit number of drones\n",
    "n_drones = 4\n",
    "# assume a 100 km/hr drone speed (about 38 mph) for each drone\n",
    "speed = 100\n",
    "# average time to prep, dispatch, clean, perform maintenance etc. a drone is 25 min, converted to hours\n",
    "maintenance = 25/60\n",
    "# calculate average service time from i to j (converting to meters to kilometers)\n",
    "distance = df['Harversine Distance (meters)']\n",
    "service_time = maintenance + 2*0.001*distance/speed\n",
    "service_time = service_time.rename('service_time')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8781d6e4-b32f-4f96-a4af-cd4f6a7a99ee",
   "metadata": {},
   "source": [
    "Next, we use the MLE for the arrival rate of an OHCA to each location $i$ as $\\lambda_{i} = N_i/t$, where $N_i$ is the number of incidents over time $t$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7ea9f343-6334-4bad-b3f2-10a29d0090b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1700 PLEASURE HOUSE RD    0.066672\n",
       "1800 DOLINA DR            0.033336\n",
       "3600 BONNEY RD            0.033336\n",
       "3400 CHAMPLAIN LA         0.033336\n",
       "2400 LAUREL COVE DR       0.033336\n",
       "                            ...   \n",
       "4400 SANIBEL CI           0.033336\n",
       "1600 OLD DONATION PKWY    0.033336\n",
       "600 FLEET DR              0.033336\n",
       "200 ARAGONA BL            0.033336\n",
       "1200 GENERAL ST           0.033336\n",
       "Name: Incident_Location, Length: 68, dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "incident_count = ohca_df.Incident_Location.value_counts()\n",
    "# arrival rate is in hours to match drone speed units\n",
    "arrival_rate = round(incident_count/(24*n_days),6) \n",
    "\n",
    "I = df.index.get_level_values(0).unique()\n",
    "J = df.index.get_level_values(1).unique()\n",
    "\n",
    "# make list of each OHCA / base pair\n",
    "base_to_ohca = df.index.tolist()\n",
    "\n",
    "# print DAILY arrival rates for each location \n",
    "24*arrival_rate.sort_values(ascending = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12edffa7-8a19-458a-b6d3-c176f7ed28e6",
   "metadata": {},
   "source": [
    "#### Base Formualtion (BF) Variables and Constraints\n",
    "\n",
    "We start to build the MILP model by introducing binary variables $x_{j}$ and $y_{i,j}$ which will be 1 if location $j$ is chosen and services OHCA $i$. We start with constraints typical to facility location problem. \n",
    "\\begin{align}\n",
    "&y_{i,j} \\le x_j \\\\\n",
    "&\\sum_{j \\in J} x_j \\le p \\\\\n",
    "&\\sum_{i \\in I_{j}} y_{i,j} = 1 \\\\\n",
    "\\end{align}\n",
    "\n",
    "Next, we require each queue meet stability requirements:\n",
    "\\begin{align}\n",
    "\\sum_{i \\in I_{j}}\\lambda_{i}s_{i,j}y_{i,j} \\le 1\n",
    "\\end{align}\n",
    "\n",
    "Below are these constaints written in gurobipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e5b07a57-a0f6-414e-8d91-0b180939db79",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gp.Model('drone_network')\n",
    "\n",
    "x = m.addVars(J, vtype=GRB.BINARY, name='x')            #selection\n",
    "y = m.addVars(base_to_ohca, vtype=GRB.BINARY, name='y') #assignment\n",
    "\n",
    "m.addConstrs((y[(i,j)] <= x[j] for i,j in base_to_ohca), name='select')\n",
    "m.addConstr(x.sum() <= n_drones, name = 'drone_limit')\n",
    "m.addConstrs((y.sum(i,'*') == 1 for i in I), name = 'ohca');\n",
    "m.addConstrs(((gp.quicksum(arrival_rate[i]*service_time[i,j]*y[(i,j)] for i in I) <= 1) for j in J), name='stability');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01193fab-c732-4f14-b0ef-445900d60e2f",
   "metadata": {},
   "source": [
    "#### Reformulated MILP: Auxilary Variables, Constraints, and Objective Function\n",
    "The objective for the original formulation is to minimize the average expected response time of the drone network. \n",
    "\\begin{align}\n",
    "\\text{min}\\space \\frac{1}{\\sum_{k \\in I}\\lambda_k} \\sum_{i \\in I} \\sum_{j \\in J_{i}} \\left(\\frac{\\lambda_iy_{i,j} \\sum_{k \\in I}\\lambda_ky_{k,j}s^{2}_{k,j} }{2[1-\\sum_{k \\in I}\\lambda_ky_{k,j}s_{k,j}]} + \\frac{d_{i,j}\\lambda_iy_{i,j}}{v_j}\\right)\n",
    "\\end{align}\n",
    "\n",
    "The objective function contains factional and bilinear terms but can be reformualted in such a way that these terms are removed and linearized, creating a MILP. To do this, Linking variables need to be introduced to create the reformulation.\n",
    "\n",
    "\\begin{align}\n",
    "&V_j = \\sum_{i \\in I_{j}} \\lambda_is^{2}_{i,j}y_{i,j} + \\sum_{i \\in I_{j}} \\lambda_is_{i,j}y_{i,j} \\\\\n",
    "&\\nu_{i,j} \\ge 0 \\\\\n",
    "&\\nu_{i,j} \\ge \\overline{V}(y_{i,j}-1) + V_j \\\\\n",
    "&\\nu_{i,j} \\le \\overline{V}y_{i,j} \\\\\n",
    "&\\nu_{i,j} \\le V_j \\\\\n",
    "\\end{align}\n",
    "\n",
    "The objective function now is:\n",
    "\\begin{align}\n",
    "\\text{min}\\space \\sum_{i \\in I} \\sum_{j \\in J_{i}} \\left(\\frac{\\nu_{i,j}}{2} + \\frac{d_{i,j}y_{i,j}}{v_j}\\right)\\frac{\\lambda_i}{\\sum_{k \\in I}\\lambda_k}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "59906f9a-a02e-46b9-a2d6-4e3312c1be7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxilary variables and constraints for R-MILP\n",
    "Vbar = 2*service_time.sum()\n",
    "V = m.addVars(J, ub = Vbar, name = 'V')\n",
    "nu = m.addVars(base_to_ohca, name='nu')\n",
    "\n",
    "m.addConstrs((V[j] == gp.quicksum(arrival_rate[i]*service_time[i,j]**2*y[i,j] for i in I) + gp.quicksum(arrival_rate[i]*service_time[i,j]*nu[i,j] for i in I) for j in J), name = 'aux0')\n",
    "m.addConstrs((nu[i,j] <= V[j] for i, j in base_to_ohca), name = 'aux1')\n",
    "m.addConstrs((nu[i,j] <= Vbar*y[i,j] for i, j in base_to_ohca), name = 'aux2')\n",
    "m.addConstrs((nu[i,j] >= Vbar*(y[i,j]-1) + V[j] for i, j in base_to_ohca), name = 'aux3')\n",
    "\n",
    "# for ease in setting the objective there is a bit a simplification\n",
    "L = arrival_rate.sum()\n",
    "m.setObjective((1/L)*(gp.quicksum(0.5*arrival_rate[i]*nu[(i,j)] + (arrival_rate[i]*distance[i,j]/speed)*y[(i,j)] for i, j in base_to_ohca)));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f91b95e2-c509-4b37-990f-e884d8b3ae23",
   "metadata": {},
   "source": [
    "Run the optimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "80838528-6c2e-4cd6-a5ff-0fed6b205829",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (win64)\n",
      "\n",
      "CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "\n",
      "Optimize a model with 1987 rows, 966 columns and 6202 nonzeros\n",
      "Model fingerprint: 0xdb44b4e9\n",
      "Variable types: 483 continuous, 483 integer (483 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [2e-04, 6e+02]\n",
      "  Objective range  [7e-03, 5e+00]\n",
      "  Bounds range     [1e+00, 6e+02]\n",
      "  RHS range        [1e+00, 6e+02]\n",
      "Presolve removed 7 rows and 0 columns\n",
      "Presolve time: 0.06s\n",
      "Presolved: 1980 rows, 966 columns, 6202 nonzeros\n",
      "Variable types: 483 continuous, 483 integer (483 binary)\n",
      "\n",
      "Root relaxation: objective 3.732321e+01, 202 iterations, 0.01 seconds (0.01 work units)\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "*    0     0               0      37.3232102   37.32321  0.00%     -    0s\n",
      "\n",
      "Explored 1 nodes (202 simplex iterations) in 0.10 seconds (0.09 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 37.3232 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 3.732321020846e+01, best bound 3.732321020846e+01, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb57b836-1f97-4fe9-88c8-a0549b82d799",
   "metadata": {},
   "source": [
    "#### Visualize the Solution Network\n",
    "\n",
    "Use `getAttr` in gurobipy to quickly get the $x$ and $y$ decision variable values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f926a13c-7915-4274-bd60-d0d6b894ba43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21 MUNICIPAL CTR Fire Station           1.0\n",
       "800 VIRGINIA BEACH BLVD Fire Station    1.0\n",
       "3769 E. STRATFORD ROAD Rescue Squad     1.0\n",
       "5145 RURITAN COURT Rescue Squad         1.0\n",
       "dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# recover the bases used in the optimal solution\n",
    "sol_x = pd.Series(m.getAttr('X',x))\n",
    "bases_selected = sol_x[sol_x > 0.5]\n",
    "bases_selected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64ed7145-4574-4066-a1d4-67615f8091f2",
   "metadata": {},
   "source": [
    "Now export the assignment variables and create a dataframe with the coordinates for each final indicent/base pair"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b305b2f6-32c6-4f82-a563-14dc646766b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# recover the assignment variable\n",
    "sol_y = pd.Series(m.getAttr('X',y))\n",
    "sol_y.name = 'Assignments'\n",
    "sol_y.index.names = ['Incident_Location','Base']\n",
    "assignment = sol_y[sol_y > 0.5].to_frame()\n",
    "assignment = pd.merge(assignment.reset_index()[['Incident_Location','Base']], ohca_df[['Incident_Location','Latitude','Longitude']])\n",
    "assignment.rename(columns = {'Latitude':'Inc_Latitude', 'Longitude': 'Inc_Longitude'}, inplace = True)\n",
    "assignment = pd.merge(assignment, stations_df[['Base','Latitude','Longitude']])\n",
    "assignment.rename(columns = {'Latitude':'Base_Latitude', 'Longitude': 'Base_Longitude'}, inplace = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e8da85a-70b5-4fb7-9d2f-b26c1230537c",
   "metadata": {},
   "source": [
    "We can now visualize the solution by finding the coordinates for each assignment pair and creating a plot with the incidents and stations to visualize the final network. The light orange points are stations that were not selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ee414c95-bd81-41b7-b154-146994052c52",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABi30lEQVR4nO39eZRk6VnYCf+eu8aekXtm7VXdVS11t7paUklqwC1BS4JGA8hgFjFwRmfgswyHxfaMMWjgmAGPviNkmJnzfXjMkWUd63w2gjZIAgu0NNitBaRu9VZSd/VSXUtXZWVW7pGx3/X9/rg3sqIysyqXyqrc3t85mRFx731v3DduxH3us4tSCo1Go9FoujG2+gA0Go1Gs/3QwkGj0Wg0y9DCQaPRaDTL0MJBo9FoNMvQwkGj0Wg0y7C2+gA2g4GBAXXkyJGtPgyNRqPZUTzzzDMzSqnBldbtCuFw5MgRnn766a0+DI1Go9lRiMjrN1q3qllJRDIi8pSInBaRF0Xkd7rW/YqIvJIu/9gNxv9TEXkh3eafdS3/NyLysoh8W0Q+KyLldPkREWmJyPPp3x+tZ7IajUajuXXWojl4wCNKqbqI2MDXReQLQBZ4P/CAUsoTkaGlA0XkfuAfA28HfOCLIvJXSqmzwOPAh5VSoYj8HvBh4NfToeeUUg/e6uQ0Go1GszFW1RxUQj19aad/CvhF4KNKKS/dbmqF4W8EvqmUaiqlQuArwI+m2385XQbwTeDALc1Eo9FoNJvGmqKVRMQUkeeBKeBxpdSTwAngYRF5UkS+IiJvW2HoC8A7RaRfRHLA+4CDK2z3c8AXul4fFZHn0v0+vJ4JaTQajebWWZNDWikVAQ+mfoHPpuYiC+gFHgLeBjwmIsdUV7EmpdRLqcnocaAOnAbC7n2LyG+my/5zumgCOKSUmhWRtwKfE5H7lFLVJeM+BHwI4NChQ+ubtUaj0WhuyrryHJRSFeAJ4FFgDPhManZ6CoiBgRXG/Ael1FuUUu8E5oCznXUi8kHgh4Cf6QgVpZSnlJpNnz8DnCPRUpbu9+NKqVNKqVODgytGYmk0Go1mg6wlWmmwK5IoC7wHeBn4HPBIuvwE4AAzK4wfSh8PAT8GfDp9/SiJA/pHlFLNJe9nps+PAceB8xudoEaj0WjWz1rMSqPAp9ILtgE8ppT6vIg4wCdF5AWSSKQPKqWUiOwDPqGUel86/s9FpB8IgF9SSs2ny/8QcIHHRQQSx/UvAO8EfldEQiACfkEpNbc509VoNBrNWpDd0M/h1KlTSifBaTQazfoQkWeUUqdWWqdrK2k0Go1mGVo4aDQajWYZWjhoNBqNZhlaOGg0Go1mGVo4aDQajWYZWjhoNBqNZhm7op+DRrPj8eagdhb8CjhlKB4Ht2+rj0qzh9HCQbMnmKi0OD1WYa7h05d3OHmgzGg5u9WHleDNwcyTYOfB7YeombweeIcWEJotQ5uVNLueiUqLx89M0vIjBgouLT/i8TOTTFRaW31oCbWziWCw8iCSPNr5ZLlGs0Vo4aDZ9Zweq1DMWBRci9m6z76F/8qPj/0DRv46D587Ahf+86r7uK34FTBz1y8zc8lyjWaL0MJBs+uZa/jkXYtYQWnuS2TO/yFeu4qgoPk6PPWhrRUQTjkxJXUTNZPlGs0WoX0Oml1PX96h4YUUMzbfU/vXnIsd/q7+AAetCUadOSwjxnzm/8De/1OYhmAZBqYh63qPW/JpFI8nPgZINIaoCUEDBu5f50w1ms1DCwfNrueNIyU++9wVoljhzxykHdlUogK1MMflcJQBq0Le92CqsThGBCxTsAzBNAwsQ7BMuU542Onr6arHf3t5imLGYqDg0vBCHj8zyXvvHV6bgHD7Eudz7Sx4s4nGMHC/dkZrthQtHDS7kpl6m6mqx1TNo9IM6MlajFVauK7LPfIqw9YsSoQxf5hW7ODm+hnozSAIYayIYkUQxUSxIowVfhgTRDErFTH++3MzeGFErBTFjE0xYwOJr2PN2oPbB+47NvET0GhuDS0cNLuCph8yudBmquYxXfMIouQqXspa3D1UYKjYz0DBwXj95xMfQ9QG4J7M68yoUaaO/w7jlTbDpQwDBYe0x8gy4lRYJEIjJkzfZ7iUCJYOeddipu7d5llrNLcPLRw0O5I4jpmu+0xW20zXPKqtpDW5bQpDJZfhUobhokvGWfIVP/ozyePp34TmJYz8IYZOfoSegz/BRKXN1YU2labPvnKWvLv852EYgrPojzABONyfo+VHixoDQMML6cs7mz5vjeZOoYWDZsdQawfXtIO6RxwnvoHenM09owVGiln6Cmu4IB/9mWtCIsUFjgzkWWgFTCy0OD/doJyzGe3JYJk3D+o7eaDM42cmgURjaHghtXbIQ8f6NzpVjWbLWVU4iEgG+CrJ78cC/kwp9dvpul8BfhkIgb9SSv3LFcb/U+AfAwL8e6XU/50u7wP+FDgCXAR+stNCVEQ+DPw8SZvQX1VKfelWJqnZmYRhzGTHd1D1aPoRABnH4EBvLtEQihkca/MisnuyNkXXYqrmMVP3qLYDRkoZ+vI3NjWNlrO8995hTo9VmKl79OUdHjrWv30ysDWaDbAWzcEDHlFK1UXEBr4uIl8AssD7gQeUUp6IDC0dKCL3kwiGt5P0mf6iiPyVUuos8BvA3yqlPioiv5G+/nURuRf4AHAfsA/4GxE5oZSKbn26mu1OpelzdaHNZLXNfDNAKTAM6M87HB3MM1x06cndXnONYQgjPRnKOZvxSovxSpv51NSUW2qmShktZ7Uw0OwqVhUOKmkyXU9f2umfAn4R+KhSyku3m1ph+BuBbyqlmgAi8hXgR4GPkQiW7023+xTwBPDr6fI/Sfd7QUReIxEu31j/9DTbnbYfMl3zmai2mKn5eGEMQDFjcnQgx1Axw2DBxdpE7WCtZGyTY4MFKk2fiYU256Ya9BUchovuqqYmjWansyafg4iYwDPA3cC/VUo9KSIngIdF5CNAG/gXSqlvLRn6AvAREekHWsD7gKfTdcNKqQkApdREl+axH/hm1z7G0mVLj+lDwIcADh06tJZpaO4Qpy/N89cvTDBZTaJ/3nf/KCcP9QKJI3mu6XN1IQkzXWgGAJgmDBUyDKamokJm+7jDyjmHYsZmqtZmtu6z0AwY6cloh7NmV7OmX2Bq0nlQRMrAZ1NzkQX0Ag8BbwMeE5FjqabRGfeSiPwe8DiJ9nGaxD9xM1Yy7C6LLldKfRz4OMCpU6dWiD7XbAWnL83z8a9eoJy32NeTZaEd8P88cY73nxylr5Bhqt4mSg2EPTmb48MFRnpc+nIOhrF978ZNQxjtydKbc7hSaXFlvsV802d/OUvGNrf68DSaTWddt2dKqYqIPAE8SnJH/5lUGDwlIjEwAEwvGfMfgP8AICL/73QcwKSIjKZawyjQMUuNAQe7dnEAGF/XrDRbxl+/MEE5b9Gbc5lr+Cy0QqpewH95doyfeOshRkoZRktZBovO8jDTHUDGNrlrsMB8IzE1nZ2s019wGC5l1l1yQ6PZzqx6qyYig6nGgIhkgfcALwOfAx5Jl58AHGBmhfFD6eMh4MeAT6er/hL4YPr8g8BfdC3/gIi4InIUOA48tf6pabaCyWqbnjTeP4pjDIFD5SxZ2+R9D4zy9qP9HOzP7UjB0E1v3uGekSJ9BYfZus+rkzUqTX+rD0uj2TTW8gsdBT6V+h0M4DGl1OdFxAE+KSIvkEQifVAppURkH/AJpdT70vF/nvocAuCXOuGqwEdJTFE/D1wCfgJAKfWiiDwGnCExQf2SjlTaOQyXMiy0A3pzLoPFDADzTY+D+dwqI3cepiHsL2fpTaOaLs+1mGskUU3a1KTZ6YhaqVjMDuPUqVPq6aefXn1DzW2n2+fQk7FZaAdUGiEfeufRRaf0bkQpxVzD52q1jVIwUHAZKroY2tSk2caIyDNKqVMrrdu+HkDNjuTkoV4+9M6jFFyL8YUWBdfa9YIBQEToL7jcM1ykJ2szXfN4darGQivY6kPTaDbEzjb8arYlJw/17nphcCMs0+BgX46+fMh4pcWl2SbFjMVoOYNraVOTZuegNQeN5jaQd5NqsCM9GRp+yNnJOpPVNnG88824mr2B1hw0mtuEiDBYdCnnbK4uJDWiKs2A0XKGUlcFV41mO6I1B43mNmOnpqajg3lE4PWZJq/PNvDTUiEazXZECweN5g5RcC2ODxUY7nGptUNenawxVWuzGyIGNbsPbVbSaO4gIsJQMUM56zCx0GJyITE17StnKazQXEij2Sq05qDRbAGOZXC4P8/hgRxKwYXpBpfnmgSRNjVptgf6VkWj2UJKGZvCkMV0Pel9XW0HDJcy9N+kuZBGcyfQmoNGs8UYhjBcynB8uEDOsZiotHltqk7DW62AsUZz+9DCQaPZJriWydGBPIf6c0RKcT41NYXa1KTZArRZSaPZZmykj7VmdSYqLU6PVZhr+PTlHU4eKOvWrjdBaw4azTak08f67qECWdtkvNLm3HSDlq8LFG+EiUqLx89M0vIjBgouLT/i8TOTTFRaW31o2xYtHDSabUynj/XBvixBFPPaVJ0rlRaRLsOxLk6PVShmLIoZG6WgmLEpZixOj1W2+tC2LVo4aDQ7gHLO4cRwkf6Cw1zd55WrNeYburnQWplr+ORdi+max9VqG0jqX83pz/CGaOGg0ewQTEPYV85y91ABxzIYm29xbrpOO9CmptXoyzs0vBDLFMIoRilFwwvpyztbfWjbFi0cNJodRtYxuXuowP7eLF6QmJomFrSp6WacPFCm1g5pBxGxUlSaPrV2yMkD5a0+tG3LWnpIZ0TkKRE5LSIvisjvdK37FRF5JV3+sRuM/+fp+hdE5NMikkmX/6mIPJ/+XRSR59PlR0Sk1bXujzZprhrNrqIv73BiuEA5ZzNTS/pYLzR1c6GVGC1nee+9wxQzFpVmgGUavPfeYR2tdBPWEsrqAY8opeoiYgNfF5EvAFng/cADSilPRIaWDhSR/cCvAvcqpVppb+gPAP9RKfVTXdv9AbDQNfScUurBDc9Ko9kjWKbBgd6u5kJzTQpNi9GejO5jvYTRcpYfemAfL03UGC1nGCi4W31I25pVNQeVUE9f2umfAn4R+KhSyku3m7rBLiwgKyIWkAPGu1dKErj9k8CnNzQDjUZDzrG4a7DAvnKGph/y2lSdqwu6udBSLNPAMMDT5dJXZU1JcCJiAs8AdwP/Vin1pIicAB4WkY8AbeBfKKW+1T1OKXVFRH4fuAS0gC8rpb68ZPcPA5NKqbNdy46KyHNAFfgtpdTXVjimDwEfAjh06NBapqHR7Go6faxL2aS50HTNo9LyMYDzMw2d/JXiWqbupbEG1uSQVkpFqZnnAPB2EbmfRLD0Ag8BvwY8JkvSN0Wkl8T0dBTYB+RF5GeX7P6nuV5rmAAOKaXeDPwvwB+LSGmFY/q4UuqUUurU4ODgWqah0ewJOs2Fjg3mmav7fObZK1yabVLO2jr5C3AtAy/UEV6rsa5oJaVUBXgCeBQYAz6Tmp2eAmJgYMmQ9wAXlFLTSqkA+Azw3Z2Vqanpx4A/7XoPTyk1mz5/BjgHnFjftDQaTd61WGj57CtnsC0DwxCd/EUiHIJQaZPbKqwlWmlQRMrp8yzJBf9l4HPAI+nyE4ADzCwZfgl4SERyqVbxbuClrvXvAV5WSo0teT8zfX4MOA6c38jkNJq9znwzYKQny4FyFstIfu57PfnLtRJHva8LGt6UtfgcRoFPpRdsA3hMKfV5EXGAT4rIC4APfFAppURkH/AJpdT7Ut/EnwHPAiHwHPDxrn1/gOWO6HcCvysiIRABv6CUmruVSWpuP7qo2fakk/xVzNiLy/Z68pdjJULSC2Id0XUTZDf0rz116pR6+umnt/ow9iydombFjEXetWh4IbV2qOPItwH63CwnjhUvjlcZLrkMlTJbfThbiog8o5Q6tdI6XbJbc8t0iprZprFY1KyzfK9egLYLneSv02MVZuoefXmHh471b+i8rEc73M6apGEIlik6nHUVtHDQ3DJzDZ/enMOVSou+nEMpa5N3LWbq3lYfmoZEQNzqhblbAxkouDS8kMfPTHLyQA8T1fZ1QgBYcdvtpK0kEUtaONwMXVtJc8v05R2qraRsg2kk0cx73a692+gueZ3UJwI/jPj0ty4t65HwxCtTFDPJfeds3Sdjm9suQsq1da7DamjhsA2ZqLT44gsT/PGTr/PFFya2fUz6yQNlqu2AphciArV2oIua7TI6Ja8B2kHMfNPn5as1au0AQ4RaOyTvWBQzFt+5skDetYgVtPyI8UoLP4yZrm0fTdIxDaJY6RasN0ELh23GTuxYNVrO8vDxQVzbZL7pk3XMbWVC0Nw6nainzvN9PVmafkTWNjk/U+fcdPLnd5XD7sna7OvNUHAtJqtt2kHEdM3bFvkFrp1GLGnt4YZon8M2o1t9b/rhjnHu9hdcvvuuAe7dV1o0LWl2DycPlHn8zCSQ5El4YUQpa3N8sEBPzma24VNtB1y90ibnmFycbbCvJ0tv3sGxDPKuxZsO9HB1oc1sw2OklKGc2zqzo5uGs/phTF7X31sRrTlsMzrqe8MLmap5LLSCHZG0FMUKEbRg2KV0op6yjslM3SPrmPz02w5iWwamIRzpzzFaypBzTd44WqI/73B5vsmZ8QVE4H1vGuXtR/s5OpjHMoTLcy1em6pRT7WRO41jGohozeFmaM1hm9GdtNT0I+abPi0/ZKC49tubrQgjDKIYy9SCYTezUtTTUCnD6bEKsw2fgaLLu984TDFrc3WhxZX5Ngstn5yT3NzkHItS1uLuoSKVps/VapsL0w2KGYuRO1xiXERwdI2lm6KFwzajW33vy9tU2wFj8y3uGszzxRcmVr3g3yjk8Hb7AMJYYWmtYc9xozDZu4eKHOzLMVXzuDTb5MJMg8vzTQaLLof780nIcyYxR03V2pydrNObtxkuZbDNO2PQcC1DRyzdBG1W2mZ0q++zDZ+DvVnu21fiq2dnWGgGqzqpOz6LgmvR9CL8MObCTJ0//O9nb2vkUxTHi7V7NBpIahgd7M3x0LF+3na0l4GCy2TV45vnZ3nywiyT1TblnM09w0UGig6VZsArV2tMVtt3pOWpo3MdborWHLYhS+/G/urb48nFPogohvFNndRzDZ+BgosXxpyfqXN5vkVvzkKQRaFyO7SIIFJkbK05aJZjGsJwKctQMcNCK+DyXIuJhRZPXpyjnLU53JdntJyhL+8wVfWYqnrM1n2GSy59eYclnQA2DdcyUSpxSnfqLWmuoT+RHcBCK+Bwf56FZsCZiSpBGN/QSd3xWWRsk9m6h1IxKhZKGfu2lmuOYqU1B81NERHKOYc3Hejhe+4e4L7RElGkOD1W4etnZ7gw06Av73DXUB7XNhivtDk7VWehdXv6YncilrTfYWX0r3kH0Jd3CKKY/eUM7TDilcka1Za/YgbyyQNlau2QWjsgVoo4htmGRzln0fTD2xL5FEYxSqEd0po1k3ctjg8XefjEIG89XMaxhFcmanz17DSvXq0zkHc41J9ot5dmm5ybrtP0NzeyabE6qzYtrYgWDjuAzgXfMg2ODOSZrrV5cbzKA/t7lm3b7bMQETK2weH+HIJBpRlQbwebXtYiTO3DttYcNOvEsQwO9uX5nrsH+a67++nLOVyea/LVs9OcGa/Rk7EYKbn4Ycy5qQaXZpubdqdvp/2ktVN6ZbTPYQfQXVmzFUQcGyhgGDBV89nXm1tx+9FylpMHynzpxavUvZBYxVSaIU3f4P13LW3Yd2sEaQmC6Vqbb5yf2ZaVODXbG9MQhooZBgsuC62Ai7MNxufbfKPWpjdrc2Qwhy0Gc82AanqDM1R0sW4xskkX4LsxWjjsEJY6qZ+7NM+FmQZZx+CekWUtthfH/MB9I/zda4k9N++a3Levh+FNrmEfxYqZepvz04nNeLtW4tRsfzp+iQdzDvcMR1yabfD6XJNnLy5QyJgc7M3iWCazdZ/5ps9g0WUg72JsMIzatUwam2yu2i1o4bBDefBgGS+IODNeI2ObHO7Pr7jdaDnLj586yGtTNRZaIZYhzDaSH9VmEUSKVycTO3EnkmqnlP3QbF+yjsk9oyWODRW4Mt/i/EydlybquLYwUsqQMUwmFzzmGj7DxQzlnL3uyCbXMqg0k37SGxUwu5W19JDOiMhTInJaRF4Ukd/pWvcrIvJKuvxjNxj/z9P1L4jIp0Ukky7/30Xkiog8n/69r2vMh0XktXTfP7AZE91tiAhvO9JHb97muUsVJhfaN91+tCeLYxr4UVIdsxNHvhkVYMM4ptoOyGcsZuseXpDYhHdC2Q/N9sdOfW3fe2KI77qrj3LOSctv1Jlv+jS9kLH55HWtvb7Ipo5TWveTXs5aNAcPeEQpVRcRG/i6iHwByALvBx5QSnkiMrR0oIjsB34VuFcp1RKRx0j6Rv/HdJP/Syn1+0vG3Jtucx+wD/gbETmhlNLxZkswTYPvOtbPV89O8+SFWR4+PkjvDZzNedeinLOZrEa0g4iZelIdczOyqcNI0ZdzWGgG1NIwWhfd00GzuRiGMNKTZaQny3zT59xUnYlK0mjItQ2KGYtWEFHM2Iym5Tg6pWQuTDeotHx6sjbHBguL/jDXSkp2eKHuJ72UVYWDSppM19OXdvqngF8EPqqU8tLtpm7yHlkRCYAcML7KW74f+JN0vxdE5DXg7cA3VjvWvYhrm3zXXf189dUZvnF+lnedGFysu9/NRKXF6cvzvHy1hm0anDxQpuGHFDMWhsh1TejXawoKY8W9+0qcvlzBD2OsnsxiT4eHjvVv2lw1a2M7t+jcCDeaz6kjfbT8kPPTiV9iaiFJoCtmLBZaAXEU8/xYBUG4NNfEMKDaSm5epqoe7713eLGHdBIBZW/tRLcZa3L1i4gpIs8DU8DjSqkngRPAwyLypIh8RUTetnScUuoK8PvAJWACWFBKfblrk18WkW+LyCdFpDddth+43LXNWLps6TF9SESeFpGnp6en1zKNXUvBtfmuY/2EUczfvzazaNbp0Km3FESKQ3154ljx1Ven+c5YhbxrUfdCKs0ApdSGTEFhFDPak+W77hogY5tUdE+HLWMn9gO5GavNJ+tY3Le/h/feO8zJQ2WKGZtqM+TM+AJ/9cIElZbPhdk6MTFZ2yJrm0ynAuT0WAWz00860GalpaxJOCilIqXUg8AB4O0icj+JRtALPAT8GvCYLPEGpRf89wNHSUxEeRH52XT1vwPuAh4kERx/0Bm20iGscEwfV0qdUkqdGhwcXMs0djW9eYd3HO2j4Ud84/wsUZcNtbvekmkk27qWwWS1zULTp5SxiZSi4UUbMgWFscIyhZ6szfe9YYifeegIj94/qgXDFtDdD2Sm7lFthzT9kC+fucql2SZj800mFlpMVdvM1j0qzaQPQ8MLaQdJLa47UddorXR/d68utGkHMaB48sLsdcdpmwZ3DRZ49xuGeOjufg70Zqk0AiYX2rx6tcZcNSnwN9vwuFxpkrFN5ho+E5UWT12Y5TPPju2Irot3knVFKymlKiLyBPAoyR39Z1Kz01MiEgMDQPdt/HuAC0qpaQAR+Qzw3cB/UkpNdjYSkX8PfD59OQYc7NrHAVY3RWmA4Z4sbz5Y5tlLFb51cZ53HOtDRBbrLSUlwAOytkFPzmG+6XOl0uZwfw7LEK5WW+Qca12mIKUUUaywTYOGF5DVdtstpXOuIe1ZgMJwLWYbPu0wIoqT86XWcP03jCT/wBRJHg3BWOm5yLVtu7bfjJpInfnEKtl/O4wIo5izk3XOjFexLSFrm2Rtk4yTPO4vZ9lfzjK50Ob12SaTNY/5dkA9iIjimJxtcWZ8gYGCw+NnJgkjRcG1bmvtsZ3IqsJBRAaBIBUMWZIL/u+R+CEeAZ4QkROAA8wsGX4JeEhEckALeDfwdLrfUaXURLrdjwIvpM//EvhjEfk/SbSN48BTG5/i3uLwQJ5WGPHSeI3nL1d486He63pE1L0QL4iJ45ihksv9+0s0/Qg/ihER/sHdA+v6YXQuNCJJpmlv7jbYbb05qJ0FvwJOGYrHwe3b/PfZBXSf606ntVo7YLDocmK4uLhdR6hHncc4KbXSeR13L0+f+2F8bf0arDCd5k83FCTdgqbrdbcA6p5PJz+n0kwExnCPS9uPaQUR1da1XAXLFDK2yYmREnONgLcd7uX05QUW2j4tL8Yg4vmxCvt7MrxxtIfj/pfJXvr/8aDzJIE9wtnow4w+/Eube2J2IGvRHEaBT4mISWKGekwp9XkRcYBPisgLgA98UCmlRGQf8Aml1PuUUk+KyJ8BzwIh8Bzw8XS/HxORB0lMRheBfwKglHoxjWo6k475JR2ptD7eMFKi5UdcnGmSdczrekT05mzOTddp+xFvPtyHbRq86UCR/6Evx8tXq9jrqE45UWnxrYtzvHK1xoHeLCM9GQ71Lc/YviW8OZh5Euw8uP0QNZPXA+/QAmIFlrbzbHjhioEBIomt/VYSnZYKj0gl+QI3EzqBitelvWRtk++MLVBwrcX2pE0/4l0nBlEqyYUoZJJZBGGMH8V4YYQXRBgCJ0YKvDpZZ185S8mzEEns1oYIL1+t4c++SF/4x8xjcq9lkQ+vcN/Yr8GFMhz9mVv4dHY+otZyhrY5p06dUk8//fRWH8a2QinFk+dnmVjweMuhMo5lLEZ8WKYwWHAZLmWYrSfO43tGitTaIXMNn3tGijdtuDJRafHEK1P8/blZShmL3rxDX85mqubzgbcf5MhAYfMmMvMkxB5YXUl+YQMMNxEQmmVsp2ilmx2Luk6oJPkyK2kvVystXpyoMt/wKWQsjg8V6Ss4q2ovsVKEsSIII/xI4YURQRijUEzXPf7bmUn8q1/BCyErbX573ye4J38lGZw7DP/w4u39cLYBIvKMUurUSut0hvQupZMk97XXZnjucoWHjvXz6P2ji+vPTddpehHFTGKPnqy2OdCbS8oSNPzFEL+ldKJHLszUGSw6tLyYMxNV3nqoTM41eWmiurnCwa8kGkM3Zg682c17j13Gjbqz3Wlu1JXw3W8cYqQnm1QNVokGoVAYImAoUIIhECvBMhQH+/Mc6MsRq+SCr4A41TyCKCaMY8JILQqDMI4Jw+R5pBRRpNKxyXH5IWRMi5MHe1mYPcdzwXFe8w7yrea9mKZQMuv01MfZvBoCOxMtHHYx3UlyTy1JkttfzvLNc7O8dLXKTM3jhSs2j94/QiEVFoNFd0WHYid6JIwVpYxNOwiIIsWl2RZv3FdivnnzDNV139U65cSU1K05RM1kuWZb0/mumIZwZb6FUlD3Aj7/7Qm+ewPFH0WSP0Nk8dEQsEwDx5JFc5FIkjBnCAjpY7ptpBQLrYBaK+ToYJ7WlQucap7hq9WTmJKoIpNBP5Pm/WSnapSyNj1ZezFZbi+hhcMux7VNvvuuAb7y6jTfOD/LO48PUshYzDd8XhxfIIxi+gsOMzWfv3juCv/wwf2EkaLaCulZwbncHQ1ztdrCMgxsy6DuBwRRzEjPjYv6bai/dfF4YlqCRGOImhA0YOD+W/5slh7bdjHF7BY635UgismlJeQLGZP5ps9oObN4MTdEEOPaxV64/uLfueBvNPopjhXVdsB8M6DeThzXPTmbgmsy+6afZ/jF32U2uoQrEXdnxvClyMIDv8sCwuSCx+SCR9YxrhMUe+H7ogvw7wHyrpUkycUx3ziXJMmdHqsw2pOhN+8SRorenEMUK54fm8exDGYa3or76kSP7C/naPkxDT8kDCIMEbww5uSB8g2P4/RYhbxr4oUxcazW1pnO7Ut8C4abmJI6voZNdEbvtsSx7ULnu+JaJv2FpOWnbRoc7s8zUHDpL7j05h16cjaljE3Btcg5FlnHJGObOJaBZRoYGwyLrXshl+eanJmocnmuhRdGDJVcTowUODaQp+FFxPvex+Hv/i0K2TxNlYHcYZyH/h2D9/2P3D1U4J6RYnrDkwiKV6/W+cZrM3zm2TFqrdV7uu9ktOawR+jNOzx0tJ+/PzfLN87PMl1rM1zKEit4fS4gY5kYhnBxtsW77rGZqnq0/Iisc7063YmGsQw4MVzkpYkFqn7I/Qd6Vo0Pn2v4SFqqo+BaWGYiuGbqKwuiRdw+cG+f87lj/rBNg1gpXVF2k1hr5NRm0g4iKs2ASssnCBWGAeWcTW/Oua6szHilRcOLONCbJbf/Z8g0v49a3SN46BB2lwnJsQwGiy6DxaTh0EIr4K9enyOKFQvtEEToTUOGd9v3RQuHO8xWqqNDpQxvOVTmmdcrzDUC8o5FOefQn3OZrnuYRpJpGscKkaS96AHn+tDUTuOhv315kqoXcvdwkXceH+SBg703LBveIWubXJ5rMtKTXSxyttXF+bww4uJsA9c0iZSinE3yA9YktDQ3pbtJ1Uzdoy/v8NCx/k3/vodRctGebwa0/CTqvZCxGCklGsnSUtxzDZ/Zus9A0Vn0wWUdE4Wi4UeUb+Bf6AgKxzK4e6hAK4gX+1Dvxu+LFg53kA3Z3DeZQ/15WkFMpenx0kSNe/eVONifZbLWolIPePvRfsYXWgwUXSrNgJFSvKzb1mg5y7tODNH0Iy7M1CFNgnNvkiMRRjEDBZdz03VMI4k6uRN3kisRxyq9mPg0vAhBCKKYoVKGnLM9hNZu4XZFTimlqLZDKk2fWjtEKcjYBiM9SV+HG4ViN7yQ8UorFR7X/GN51wSVrO8kD96IvrxDy4/oyV7zye3G74sWDneQ7ro3s3UP1zbJu+YdV0fvGSnS8kNgnrmGR9axuHuogB8mjsPJqkd/3kUpmGv6DBVXdjIrpWiHMRnLwLGMm9qFJxba9OYd/tFbDvDKZO223kneiLoXMt/wWWgFKJXcCQ6XXH7w/hGeeGWaTs6Prii7fWn5EfNNn0ozIIrVYhZ1b85ZZgJdih/GvD7bxDYNDvXlrvu+Zhc12dXzbbfCXLYVaOFwB+lEb4RxkvJf85LIiYmFNtV2QNG1NqUezVo4ebBMO4iZWGjz5kM9HO7P843zs0k7URHGKy36Cw5zDZ/BwsphrX4Yo5I6aGRuEupXSX/MwyWXoVKGo4ObmAexCl6Y2KDnmze3Qb/3XuO2mz80G8MPYyqt5DvkBTEiUMrYlPP2mn8zcay4NNdAoTjcn8dcYmpybQPTSPxhq3GnzGVbjRYOd5DuOjEHenO0g4ipapuMbfD6TBPTEErZxA9QWKEnw2aSJMn18rXXZnj+8gIZ2+LBA2Xm0t68oqCYtUFBtR1ep0J38KOYWClEBNdeWY0PopgrlRZZx9zU1qQ3I+oyGzW91W3QsH0SxzQJHdNfpXUt/DTnmuzvzdKTtZdd3FdjbL5Fy485PJBbsamPIQaubdAM1lapZy98X7RwuIMsVUeDKLHnv+eNQxSyNgvNILnLbQRYplDO2ZSzq6vLG8U0jTQHYmoxSe6+/T08eX6Wajtgptamr+AwW/dWFg5hkpnqWnJDzWEsTX460Ju97VrRUrORaxsM97iUs85iO0jN9mbpObQtYajkUs5tPBFtqtZmoRUw3ONSyqxcGFIk0X6b/uqaw15BC4c7yGrqaCljsz9N2FloBczWfWZqSQvEctamlLU3vZWhYxnXkuTOzfLw8QEO9OZ45WqVyapH1jExMGgH0bL3bocRQRTjWOaKmsNcw6feDhktZ25bC8aVzEa9eYfenE3O0V/vncBaw083QrUdMLngUc7ZN/SdQVKML2MbtP0kB2cl7XKvoX89d5jV1FHDEMo5h3LOIYzixYiMyaqXXqwNerLOTSMy1kvetfieu/r5ytkki/ptR3qZqXtcqTQpZSzMosFsw2f/kuMOwpgwVri2sSxSyQujxaiQTkb1ZrERs5Fme7FS+Glxk89hO4i4PNck6xjLvrtLESHpYhiEeGFEVt9YaOGwnbFMg768Q1/eIYhiKs2AhZbP1YU2Vxfa5F2Tcs6hlLGWhZuul57ctSS502ML3DtaTJsBtbBNg6xjMlLKXGfrnay2+c6VBeYabaZr3mLOhlKKsfkWIqz6o1wrSqnFdqZLzUa9OWfTBKXm9rHR8NONEEZJZJIhwqG+/BqEjZBzLCrNkIavhQNo4bBjsM1rmZrtIKKaOuuuzLcYl+Suqyd7a3ddQ6UMbz3Uy9Ovz+OYBkf7sjx7eYGZukch7S3dcSpPV9s8d6lCpBRDpex1XbQs06CZZp/eqq1/JZODNhvtLJp+ItS7w0/7C0n46e0wNyqluDzfIohijg3m1/wd7E7M3Gxtdyeif107kIyd1J4ZKmVo+dFimF+1FSLSoidr05Nbe5hfNwf7czSDiDPjVYYKNqM9Ga7Mt8k4Jr15Z1E4nLlaxTQgCsEUFktOfOviHMcGC5Sy1mL26XpZyWxUzFiMlhyKGUubjXYAmxF+ulEmFtrU22FSGmONNxCdJE4xoOHr3mKghcOOJ+uYZJ0soz3Z1OySRHpUmgGmIfTkbMpZe12OvXtGirSDkPPTTfryDlM1j8uzLfpzLvt6A0qZJLIqihV+HNNpF5VzTM5MVLl7qLhuc9KNzEa3w+SguT0szTyHWws/3QjzaWmM/oKzrpsTAUxTcEyD5hpyHfYCa+khnQG+Crjp9n+mlPrtdN2vAL9M0s7zr5RS/3KF8f8c+H+RtAP9DvA/K6XaIvJvgB8maTF6Ll1eEZEjwEvAK+kuvqmU+oVbmuUeoeBaFFyL/WVFzQtZaAbMN3zm6j62JfRk1x4ae/JgLy0/ZrzSoi9vc2G6ybnpKvv7sskdYM5moRWSc0xydvI1mlhok3Ms9vdm1+wD6ZiN5ps+YaTNRjuRlTLPbzX8dCM0/ZAraRDE6E1Kx6+EiCAqyZRuas0BWJvm4AGPKKXqImIDXxeRLwBZ4P3AA0opT0SGlg4Ukf3ArwL3KqVaaW/oDwD/EXgc+LBSKhSR3wM+DPx6OvScUurBW5zbnkVEKGXs60JjK83lobE9q/x43360l6+djfAWIrK2wcXZFgenG+wvZ3nDSIm/eG4c1zLI2MJMPXGSf/+9wyvmRHQTxYpK01+MVBFJBFtvj0Mpe+eyxDUb53aGn26E7tIYBzeQU5Nsrcg5JvONAJUmd+5lVj2LKik4U09f2umfAn4R+KhSyku3m7rJe2RFJABywHi6/Ze7tvkm8OMbmYDm5iwNje1knV4LjTUp55ImJktNN5NVj1or4MUrVWrtpMPbc5fmuGekSE/OZrjkks1YzDR8vDDme+4e4OSh3hWPo2M2mm8EVNvBbY1U0dweloafdoT6VvuCOqUxYqU4OpDfUOSeCCgFecdiuubhR/Ge7P7WzZpEvIiYwDPA3cC/VUo9KSIngIdF5CNAG/gXSqlvdY9TSl0Rkd8HLgEt4MtLhEKHnwP+tOv1URF5DqgCv6WU+toKx/Qh4EMAhw4dWss09jyWadCfNlnp1KZfaPlMVNpMVK6FxvZkbaaq7cUKsu841s/TF2Y5M1Gl4YU8e3GOu4YKSYb10X7efqyf2brP0cHlNWuWmo3WUyhNs/XcyfDTjXKlcvPSGGtBUt0h55oEkcIPtXBYk3BQSkXAgyJSBj4rIvenY3uBh4C3AY+JyDHVKW0JiEgvienpKFAB/ouI/KxS6j91bfObJD6L/5wumgAOKaVmReStwOdE5D6lVHXJMX0c+DjAqVOnFJp10d3EpB1Ei07sK/Mtxistnrs0j2kkCXKxUrzlcB9NP+TlqzW+dm6GjGPgRzEZx1x0AHbqQd3QbFROcjL2urq+E2j6IfPNYDHwwDJvb/jpeujuiWIawlDR5b79PTcsjbFWFIlwABZroO1l1mUcTB3GTwCPAmPAZ1Jh8JSIxMAAMN015D3ABaXUNICIfAb4buA/pa8/CPwQ8O6OUEnNVB1T1TMicg44ATy90Ulqbk4nNHa4lKHphyy0AqZqHqWMxbmpOk0/pC/v8sCBHiZrHq9fucDfXvxPxL5J0HoZ980fYHj0pxLfxm0wG+2Ffr3bgZXCT28lLPp20N0TJedYXJptcLXS5g0jpVvab7dZSSRxsu911hKtNAgEqWDIklzwf4/ED/EI8ERqYnKAmSXDLwEPiUiOxKz0btKLvIg8SuKAfpdSqrnk/eaUUpGIHAOOA+dvbZqatZJzkh/dvaNFKs2AOFZ4YVJZ1Qsi3pw5w9MTl/mWd4QBc56sdwX7hd/kVQXhyPs23Wy0HRok7Wa2Q/jpeuj0RMk6JuOVdtq1b3N6oigUtmlgm7KYX7OXWYvmMAp8KvU7GMBjSqnPi4gDfFJEXiAJR/2gUkqJyD7gE0qp96W+iT8DniUxHT1HagoC/pAkPPbx9I6kE7L6TuB3RSQEIuAXlFJzmzZjzXWsdFc+VMpw374eHj8zSd4xOdyfY7re5vXpJrnKV+iNe/lW8AauBP18o3ovl4I5Bp59jP5H3kNfzqEdhMzUY8y0MXzHntt94ynLnrDidl9/bfraJgrd33kd3EjjUipph7lS46PVIti2mk5PFEOE3pxN1k56n99Ki86JSouvvjrF1arHG0eL1NrBmkt372bWEq30beDNKyz3gZ9dYfk48L6u178N/PYK2919g/f7c+DPVzsuzXKUUsQqsfnHShHFikgp4rj7OYvLJhZafP3sDFnbxLVNXpus8/TFOe7f10PWsTBE+M74QlpCHHqzDhdqea6qYSIERcyz3n1MxlMUvDbZM1NYhmCZknaHM3EsA9sysA3Btgwc08BKk41sy8AykiYrppEIB8sQDOOaQLk406Scs/EaMRnbxEB2Zb/eDeHNQe0s+BVwylA8Dm4fsLLG9dffmeDUkT5c29gW4acbobsnSudGodYONtyis/M5hVFMOWvjBTEvjtewDAN2WWe39bIzvhF7AKWWX8CjOL2wq2uPybKu9Z1xsUKtwS0vck14fOvCHHEc0/AVk9U2QaRo+iHfOD/H/ft7cG2DNx/sRSlFK4yotULmXYOesMabnHMcz1xi1KnQY9XJZArkDpYwxcAAglgRRDExibqZvjuxgiBSeEGEUiHVts/luRYtP6SUtblrsMBgMYNtCa5lUspY+GFMT9ah2kqyvpt+SN61qLYDLEMS4SKyqKnA5vkptq2/w5uDmSfBzoPbD1EzeT3wDnD7rmtJW20FNPyQhhfy5PlZvv++kS0PP90om92is/M5xSqJrCtmbHqzNmen6vhhvKf7gOx54bAZP/7rLtyrXMDj65Zdu1Cv9cJuphdDI70YOmJgGCxeIA1DMIBYQRgnzXhCFRNFiiBWSRayJOObfkQ55xIphWUaZCwT1xYa7ZB3nRjAD2Pmmj4zdY+4BW7BIDzwTgbPf43DzgRvK5yhrVwKZoxz788SZBz6ckljnbxrkrGSsMBWEOEFEVGsCONkoraZmAJemajSk7MYKmWotwNeGK9y8oBBOW/T8mN6sjbPXa6QsdtkbBMviPHDmLcc6eX51yuLn4dpJFqHa5vMNTyePD9HMZNkjI9XWpybqvO99wyxvze7ZIxxQ7v6tvZ31M4mgsHKJ687j7Wz4L5j0fwC15yr+8pZWkHEkYH8VhzxprDZLTq7zVSd5M2eXCocIi0c9iydH3/BNSllbCoNn/96epx3nhhksOguM8N07uzVBi7sSy/glmHgmLJs+eLj4nOW3RVDomn4UXKh9MMYL4xpBtHisu5j6hQVy9omTtbATc09bxwtEkRqUT33goir1TaubTJV86i2AtJrOft7cyy0PMKekxy8Cwoz/5W77AnmreN4R/8xfXc9CiS2bCOCehsaEtGXdzjSk8MykoZBrSCi5Ud4YcRTF+bIZyxGStlFs0CtHeCmDYjCWBFGMW86UObbYxVmGz6jPSbHBgsUMzZhHC+ehyCMk/PRCnj29Tn8MMbwSKOkhChWfP21Gb77roEVz4+ZmsMMkUVt5Gtnp4nTD1JtN3+HX0k0hm7MHHizwPXml5FSBsMQau3gjrVqvZ1sZovO7s+pQxwrsraJH8aJV3SPsqeFQ0eldC2T8YUWkKiWX3l1+rqLyEoXdtfqvphfu4AbXSYOY/GRDYUBdgRAK4jwwmuC4GYCwLUMShkbxzJw0tc3CiN96+E+vvjCVWrtRAhUWwHVdshbDpWJYkXetfCjmOEeF1GKF660KLgWB0/8A/y7voe+4wNYrYCJaptYKYZLGbJeSNOLcCzBtQzm0kJoxYzFQNG9rhTyk+dn2VfOXnfn3vEndMfSl3MO9+3vWXb8cawIOtpRlDyPYsVLEwuMlDLECkZ7MogIURwzVfM4PlxYFPJRlGgysUoeoyhZHsYx7VAxXkkq3M42fFzLxDS2kb/DKSemJKtLC4iayXKWm19q7eCWzC+7lZXMVEEUc7gvR9MLN+zL2A3saeHQrXoPpqolRcVc0+cNo8VFYXA7USoJFV2qBfhhTBAtFwAZOzH/dASAmwqB9eQRdLJe/Sjm6ECeVyZr1NsBhYzFd+3rY7CUIYqS2jIDRZeCY/Hll65iGsJ9+3po+RG2CbZhMFB0aXgRsYqZrnncNVig7oVMVttEccxoOUMYKWbrPhemG2Rsg4FCUpStv+DS8iMy9rW7tsY6fpCGIbiGyVJf6l2DBVp+dN3dYNOPGCy660rgujjToOmF5F2bztdgPcd3WykeT3wMkGgMUROCBgzcD2y++WW3stLn9N57R3hlsrYY2rtX2dPCoVul7ERr1NoBw6XMppYFWCoAurWApQLAMK6ZgMo5G8c0FrWAWz2mdhAx3/SZb1zLej02lOfQQG4xrjtrm7TDCEOSzNNixuTxM1M02hFv2l8iY5uEcYyIYFtJ9yzXNlAqcTZfmmty91CBYsbi8lyT8fk2vXmb40MFal7IbN1jbL7FxEKb/eUsT1+cBzbHudhhs5yWnf2IJBpDfTvdfbt9ifO5djYxJTnlRDCk0UqwueaX3czSz6kdRJybadAMgi08qq1nTwuHzYx8iOPEBOQtMf14YUQQXu+U6AiAnGPiWNcEgGsZt9zucymdxjlzDX+xlEUxY1HKWLSCiLlG8gMo52xQMN8McCyDg31ZXMvk716bYbzS4v79JYpZGwPBEIOMnYShAvTnHcYrbYZLLlM1jyvzLQ7157h7qMBk1WO65tFIO8MdH07iyGfrPrV2UihtrNKkmgrlzbi73ay75m1/9+32gfuOrT6KXYdtGtiWUNeaw96l+8d/bqpOpeXTk7U5PVZZXN/NSgLAS53ASwWAaUgSteNYOLlr5p8kzv/2R0AsrbHfaZxTylgstALGF9rEcSIUyjmbyapHy4/ozdvs68kiAt8ZW+DsVI2DfXkO9edZaIZkHIO4GeNYFraZ2Fp6c04yPogYLmW4utBmpu4xUHAZ6clQzFiMzbc4P91gsOgyXHIpZmzaQURfwWGg6CalC1yTfMbalAiyzbpr1nffew/TEDKWSTuICKP4jvxetyN7WjjANQEwVfUY6cmQdUwWmgF/eXqc77mrn3LeWRQENxMAbu6a+edOCYCl+GG8WPDOD+PrGudkLJPZhs+56QZRrChlrbSWUsTrs01E4FBfjp5cYqe/ONPghfEFiq7Nmw/2UPMicq5JHHdKDBiLTnYjLZkxXUs+w4ZncXWhnTQCcizyrsXdQwUmFlpM1zxq7YCDfUkFzf3lLCOlDLMNj7mGzzMX53j64jzDJXfx+LZN+Khmz5BzTOqpX04Lhz1Md8LQ5fkmUZwkaf3duVkePj6Ia18vADqhoNuh7kycNvOZbwbU20k8e941GS4lHdtEEsf767NNwkhRyFgMl1wc0+BKpUW1FZJ3TQ725RZ9GtM1j9NjFaJI8ZZDJXKOxWwjYKCYYXy+vSgcuukvOMzUPWbrPgd6s7w2XU/8D4MFLDP5rA705ihlk8qvr03VGS5lGCg4aWXNDIMFlxevLFDMWEQKDJHtFT6q2TPkXYsgauGHMbltEH+wFWjhwPVRS71ZBxEwShkWWj737ru1ao+3i5YfMdf0qTR94hhsS65rzaiUotJMqqv6YUzONTnYl6GQhjWena0TxYqRnsx1se8NL+T5yxWqrYATIwWODBQYm2/hWAZWKgyttDhZN7Zp0JO1mWv4DJcyHOrLcX66wdh867qkq1LGJjdkcqXS4upCm2o74GBvDsdKNBEvjDk2WCDu6sS1bcJHNXuGvGMSxklfh72KFg5cH7VUyFyLWuqOyd8OhFFMpZX0hW53lVQu5+zrwjYXWgGT1TZeEJN1DA4P5ChlbOI4id2frSetQo/056+rnOqHMc9frjBTa3OgN8ddQ0WiWNHyI0bLGVp+8p4WsmLk1EDBpdJMnN+DxcTfMFFpM13zrhNAlmlwuD/PfMNnfKHF2aka+3qy9OadFZOStk34qGbP4NomhkDD37ulu7VwYPPrtWyUlRyxIz0Zal54XZ+ErGOwr5yhnHOuM23V2olQaPkx1ZbP1WobL4zpyzvcM1wkjBXtIKa/4CxmzXaIY8VLE1XGKy3KOYcjA3kGCi6vzzYwDOjLOVyYbWCbBn4Yrygcso5J3jWZbXgMFBwGCi5NL2KymvgflhZ3680nBd/G5puMzbeotgPu21fiv7+ctATZynOh2dvYaQRhw9+7EUt709OyhE7UUtYxmaknfZXvtAO0U8qj5UcMFFxqrYA/f2aMr5+d4fWZJnUvpL/gcHy4wN1DRfoL7qJgaHgh56brXJxpEsYK2xRem6pjiDBQcJmqtvnTb11motLi8ECOfeXssuS+89N1Lsw2yFhJLf8DvVm8MKLaCunPJ3f9LT/CtZNxS81KHQaKLkGoqLaSO679vVls0+DSXJMwWq6iO5bBscFCIgTbIbV2xHcd69/Sc6HRdMLKm3u46Y/WHFK2OmRxsYmJbTJV9WiHEbFSXJit8/6DB1Zsr9nyk7vyWjvEMoXRcob+vMOXXrxKKWtTcK3E5xAp+vIONS9csZXi1YUWr07VMUXY15PlcF8e2zQYr7QQSZzNrSBCKbBNE4humJDXydyernv05JJmMYf6cpybrnN5vsXRGxR9Gyy6i4lz7SDm/v09jPZkt4XTX7P36DT9aQcxcax2XPXazUBrDtuEuYZP3rUWL4a9OYe7h4tYRuLo7RYM7SDi0myT16bqNPyQ4R6Xe4aLDBRcRGRxXyJJqG1/PjETLbSWZ3wuNANenKiiVMxgwWUobfgSRjFzjSTvwzaNRdurY3ac0jf+sQwUHFp+RCO968o6JvvKWertkKlq+4bjMrbJ3UMFBosu840gmd8evnPTbB2dXIcgSsLY9yJaOGwTOo5YEWGkJ0NP1sYLouscsX4Yc3muydnJOtV2wFDJ5Q0jJYaK1/sPOvuCRMgUM/aKTt2WH/HSxAJtP6KUcSjnbfal2tNc00cpFh3JTS9Ky2QkY23jxl+d3pyDYcBs3b/umDrJdjfrz9uZ/7HBRMM4P93g6kIbtZbStxrNJpJzLMIoKX2zF1lVOIhIRkSeEpHTIvKiiPxO17pfEZFX0uUfu8H4f56uf0FEPi0imXR5n4g8LiJn08ferjEfFpHX0n3/wGZMdLtz8kA5tbkHxEotVtE8eaBMECU9nF+drLHQChgoOrxhpMhwKbOi2eVm++oQRDGvTtaYqfuUsjY9WZsDvTlMQ1AqKZRXyFhkbDNtKxmSc0yCWC1Wn70RhiH0510WWgFeeM2ht7+cxbUNLs81CVa5G8u7FseHCvTmbaZrHq9N1Wnr1o2aO0je3dvhrGvRHDzgEaXUSeBB4FEReUhEvg94P/CAUuo+4PeXDhSR/cCvAqeUUvcDJvCBdPVvAH+rlDoO/G36GhG5N93mPuBR4P9J+1fvalZyin/fGwYREV65WmO+4VPO2ZwYLjLak71p1uZqDvY4VlxIcxB6czYFx2Kw5FJIo4kqzYAwUgwUEk3DC2PiGPKORRDGONbq9tf+QpIv0q09GKn/IYoVl+eaq2oDRpo4d3ggRxgrXpuqM13ztBahuSNkbJNYqetucPYSa+khrYB6+tJO/xTwi8BHlVJeut3UTd4jKyIBkAPG0+XvB743ff4p4Ang19Plf5Lu94KIvAa8HfjGeia2U1gpfHWolGG27jFd94jjkHLOZqjkrqvx+80c7FcqLS7MNOjJWmRsi0LGYqSUWVyf9FMwFnMNOiaqTjLaWqrDLk2K62g4nZIZY/Mtpmoew13veyNuljin0dwu7DSzv7VHw1nX9OsSEVNEngemgMeVUk8CJ4CHReRJEfmKiLxt6Til1BUSjeISMAEsKKW+nK4eVkpNpNtNAEPp8v3A5a7djKXLlh7Th0TkaRF5enp6ei3T2HYsDV9teCGfffYKf3d2msmqR8G1OD5c4GBfbl2C4WZM1dqcn66TsU0KblJu+2BfbtHhXWsHtIP4ugTAph9hmYlzO4gU9hovygOFpKDeXMO/bnlv3qE3bzNV9ai211YWuZM4d6A3SzuIODuVaFMaze2iUydtrybCrelXrpSKlFIPAgeAt4vI/SQaQS/wEPBrwGOyJNYy9SO8HzgK7APyIvKzq7zdSjaLZXYEpdTHlVKnlFKnBgcH1zKNbUd3TaemF1Fth4RxzIXZBncN5Tncn19Xc5rVqLYDzk838MMktNU0DIZLmeveY6buY5mSlPBOqXshecdabJVqrzGsL+uY5NKkuKWmoH09WTK2wdhca1023d68w/GhIlnbZGy+xeuzjRXzJ243E5UWX3xhgj9+8nW++MIEE5XWHT8Gze0lKaAptPx4T5oy16WXK6UqJOafR0nu6D+jEp4CYmBpg973ABeUUtNKqQD4DPDd6bpJERkFSB87Zqkx4GDXPg5wzRS1q+iEnAKEcYxtGBwZSHIMcs7mpqC0g4jzU3UqzYCRsotCkXfN68patIOIejtM/QWJAPDDpA1nzjUJ4uQivJ6mQwOF65PiOhiGcLAvh0JxaQ3+h26WJs69OllfMUz3drFU42ullWO1gNhd2GZSJiZIy/TvNdYSrTQoIuX0eZbkgv8y8DngkXT5CcABZpYMvwQ8JCK5VKt4N/BSuu4vgQ+mzz8I/EXX8g+IiCsiR4HjwFMbmdx2pzvktCdrM9KTIYrVptcRCqOYC9MNpmptRkouplyrktrNdM1DJCmV0aGZqtQF1yKIkgv4zXIcllLKWItJcUvJ2CYHyjlafsTVm+Q/3IjBosvdQwVsU7g022Qsrah7u+lofJZhUG8ndaCKGWuxD4hmd2CZSR+WMFJ7MtdhLbeAo8B/F5FvA98i8Tl8HvgkcExEXgD+BPigUkqJyD4R+WuA1DfxZ8CzwHfS9/t4ut+PAu8VkbPAe9PXKKVeBB4DzgBfBH5JKbUrPULdIacKVgw5vVWUSu7Mxxda9OaSLOQoVuwvZ69z6AZRzEIroDfvXBcJ1fCjxc51HfPNejQHEUkyrP1oUdB005Ozk3LfNX9Dd/93MnGuY0r64gsTvDpZ5UqlxVzTpx1E5F1rmW9Fs/PJuUlb3L0Yziq7wZZ26tQp9fTTT2/1YWyIzeh6drN9QnLXP1BwGS1naHgRPVmbg33Xaw2T1TZTVY8TI4XrnN9nJ2tYpsHRgTxTtTaTCx737iutq6xFHCteulql6Noc6s8tW6+U4tx0HS+MuXuosGHne8MLGZtPfBidjnNLS46sl85neT7tT3HPcCkpWeIFRHHiOyllLYquRT5j8ej9o7f0fprtxcWZBmcmqjx4sLyYILqbEJFnlFKnVlqnayttMZtd06ljDy9mLFzL4OJMk/mmzw89MIofxlimLPuSx3GS9FbKWtddmMMoph3EDOdsJiot/vsrU4xXWlypNNclxDpJcdNpb4mlIagiif/htak6l+ea3DVY2NBFvZM4N76k49x8w9+QAO7+LKutENMQzk7VOdyXpdoOMQ1Fte2BKBaaAT/5toOr7lOzs7AtA4E9qTnoQPFdRsce7lomc/UABEZKLmen6vihWsyChmtmkk98/TxfOzu9LOqnmWYk11oBj5+ZpOGFG3bA9uXTpLjGyk17XMvkQG+Olh8zvrB+/0OHpYlz3zw3y2efvUJzA8feHU1W8wLKWYdYxbw+1+TkgRJF12am7jNczPDAwZ7remNodgdOJ9dhD2bna81hl9HpaidAjKKUsShlbS7NNhm411nMgu7cFRdcE9sw8MKIv3ttlpxjMVrOMlFp8cSrU4zNtwjCiNGeHL05BwUbat3pWNeS4oaKK5f96MnaDBQT/0PBsRb7WW+ETuLc0xfnCOOYph+Tc9W6jr3zWTa8kCBU1NtJSO/EQhsQ7hkp8uChMj9w3wjnZxpcqbTIOuam5aRotp65hse3L1d45vU5Ls6UefBg754pH681h11GJwJqoRVgGkJfzmWm5tFfcBguXstGXoy4MQ0ipRjtyS5G3HQEx0IzYLiYYb4Z8upkjVglRfWADTlgBwoucbw8Ka6bkVKGrGNyeb55y7WUFntX9+WIlcLgWtvRtRx757N0LIP9vRmm6h4NP6S/YDNZbTNRaXPyQDkxi6WRX2PzrT0ZE78bmai0+NqrM4SxouBa1NvhngpZ1sJhl3HyQJnpmsd4pUXeMZlretTaId93z9B1xfK6S4QXXWuxU9tcap8vuCaOaZJ1TAYLLoYBF2fri+M30rrzZklxHUSS+kuGCJfnmsS3GJral3cQuK7B0VqPvRNN1g4i7h4qcrAvS6UZMFxMWpq+cV+RXKqJOZbB/nKWphcxXdP9rncDp8cqiy14YwVZx9pTIctaOOwy+goO94wUKWYt6u2AWMEPn9zH0cHC9duld8WuZdKf9oHoXDTnGj6FjM2B3iyljM2RgRxxnORB3KjK61q5UVJcN45lcLAvSzuIGV+4tbu0tVSovRHdBQwXWgH37+vhp99+iO+/f4SfOnWQg315Ls81F2vvlHNJWfKpmrdi2K5mZzHXSCoWH+rLJhWJo3hPhSxrn8MuoBNuOVNPooHuHizw4289yOuzTXKOuWL3tZv1zT49VqHhhYv2+b68y/HhApPVNjN1j768w0PH+jdkey1lLBZaPp99bg7HMm4YPVTMJMUGp6oe9XaVi7ONDYX7di7wnc9nvce+NJqs4YVcmGnw+lyTg305Ls42uDjb4K7BAo5lsK+cpeGHXJ5rcfdQQXey28F0bqAKroUhQhirDWnMOxWd57DD6Q63bAURMzUP1zI4ebBMMWtzfKh4w+qlN8qx6N5nt+DYjF7OE5UWn33uCmEUc3QgTxirG+5bKcVTF+Z44pUpjvTn6U1/rJt1LBul1g54fbaZVJjtyXB+toFjJiU9TCPRwM5PNyjnlueTaHYO3b+DSjMgjGNca3f1NL9ZnoM2K+1wOo5lEaHpR+wrZ8nYJqfHFpZlQS9ltJzl0ftH+R/fcZhH7x9d/MKv1g/iVo93uOhS7op8upEdV0SYqSfVadthjKyy/Z2imEku+u0gYqLa5kA5ixfGizWi8q7FUMml0gyoNPeGCWI30v07qHsBlim7SjCshjYr7XA64ZaQ1ERybYP5ZkAQxZRzG1d/Nzs5r0PneHu6jq3TJ2IlFlrB4t13JzHuZtvfKXqyNnHal8KQgNGeDOOVNuMLbfaXswwVXWrtkCuVFjnH0r0ndiid38Bcw2c81bS7l+9m9Dd2h9OxixoiFDIWM3UfP4xW9DNsB7qLDXa4mR23L+8Qxeq6KrXbxe7bm3fYV04qwzb9iP6CzVzdZ6bupVnfyQXk8vz6qs5qtg8d01IYxfRk7T0VzqqFww6nOxpnruGz0PSxTYO3HOpdffAWsN7ooVuJNroT9BdchnsSE1KsoJS1mKgk3epcy9ThrTucjtm2L++ScxIf3FabNe8UWjjscDp2UYDLc00Gii4/fHLftlV71+vPuJ3+j81iqJhZrAprmULWMbg0m4S46vDWnU0nHyhjmwyXMtimsWfCWbXPYRcwWHR5w2iJ+/b3cPdg4bpkt+3Iev0Zt8v/sZmM9GSIVVLAsD9vE8ZqMcRVh7fuXDpm0E5YN2wfs+btRmsOu4DxSpsoVhzszW17wbCb2VfOUs7ZzDaCpK2qUlyaayDAwd4cfhgzvgds1buJ7W7WvJ1o4bDDmW8kTXKGSq6uCroNONCbpSdrU2kGFDM27SAJcc055mJ460LzzrU01dwaO8GsebvQZqUdjB8m5SVyblL/SLP1dKKUotmkx0PBtai1QyYW2oymPa/HKk2yzo2TEzXbi51g1rwdrKWHdEZEnhKR0yLyooj8Tte6XxGRV9LlH1th7D0i8nzXX1VE/lm67k+7ll8UkefT5UdEpNW17o82b7q7B6UUl+ebQGKyuNWOZ5rNQ0Q43Jcj55o0/BDXNpit+8w1fA72ZVFKh7dqtj9r0Rw84BGlVF1EbODrIvIFIAu8H3hAKeWJyNDSgUqpV4AHAUTEBK4An03X/VRnOxH5A2Cha+g5pdSDG5rRHmG67tH0Ig703jwLWrM1GIZwpD/PhZk6LT/CNA3GK20OD+TYnybPTdc8hkqZ1Xem0WwBqwoHldzedGo12+mfAn4R+KhSyku3m1plV+8muei/3r1QklvenwQeWd+h711afsRU1aMna9O7B6ImtiNr6f1tLgqIBq0gQgQuzTa5e6iwGN5ayFjXJfhpNNuFNd1yioiZmn2mgMeVUk8CJ4CHReRJEfmKiLxtld18APj0CssfBiaVUme7lh0VkefS/T58g2P6kIg8LSJPT09Pr2Uau4I4VozNNzENYV9Z33VuBZ2s2ZYfrdp61DINjgzkcS2TOFZEaYjrYNHFMoXLcy2iW+xZsdfptLv94ydf54svTOyJ7OU7wZqEg1IqSs08B4C3i8j9JFpHL/AQ8GvAY3IDw7eIOMCPAP9lhdU/zfVCYwI4pJR6M/C/AH8sIqUVjunjSqlTSqlTg4ODa5nGjqbzA/ijr7zGf3t5CtsQLFObk7aC7t7SQRTjWMZNs2Zt0+DoQB7XNoliRdOLGJtvsr8nq8Nbb5H1CGrN+ljX1UUpVQGeAB4FxoDPqISngBgYuMHQHwSeVUpNdi8UEQv4MeBPu97DU0rNps+fAc6RaCl7ls4PYL6RlMawDeHvz83qH8AW0cmaBag0AyYW2vhhfNNigI5lcKQ/T8Y2USTNjuaaPoNFR4e33gLPXZ7HNoUgSrSy7VC1d7ewlmilQREpp8+zwHuAl4HPkfoJROQE4AAzN9jNUu2gw3uAl5VSY0vez0yfHwOOA+fXNp3dSedO1TQMHMvkQF9O/wC2kO7igYMFl6JrMVlt4wUxtfaNL/IZ2+TYYCIggihmruETK0XWMRmrNPHD+E5NYUfT8pNaVRdmGjx5fo5nXp/jiVenePri7KLg3gvlLW43a/GEjQKfSi/YBvCYUurzqanokyLyAuADH1RKKRHZB3xCKfU+ABHJAe8F/skK+17JD/FO4HdFJAQi4BeUUnMbmdxuoVPm2hAhjhWGyLYoW71XWdpFz7EMso7JG0aKXJxpUs7ZjPZkVjT7ZeykM9/56QYLLZ+pKgyXMrSDiMvzTY4N5HVY8hKCKKbeDqmnjZ46PppaO2C+6ZOxDI705fCimOcvVzg+VGBU++NuGd0JbgfwxRcmaPnRdfVdau2ArGPy6P2jW3hke5eVopWGSxmm6x7TNQ9DkoCBG/XUqHshF6brzDUCyjmbvrxNpRky3OMyVNxdF7a1RHZ1E8eKhp8Ig3o7pB0kGpVpCMVM0pnw3HSdr52dpuWFhDGM9GTJ2AaVlk8UK/7pu0/sycS19XKzTnA6hm4HcLN+z5qt4UZZs8OlDD1Zm7H5FpfnWlSaAftW6MhXcC0OD+SJVYOZhocpQsY2mKomne92S3hrd6vNgYJLw0v6ISwtQdHyI2peQD3tjaEUiCTf95GcQ8E1AeH8dJ2/OTOJZQoLrRDbBC+MCKIYP4ooujalrKUFwyawO76Bu5xOfZfTYxVm6h59eYeHjvXrH8A2JWOb3DWYZ7bhc3WhzauTNUZ6MvTnnetMRqWMzeH+HFGsmKy1GS65QBLeenxo+1fXXQunxyq4loEhgiGyqP0+e2me77bNZaaijG3QX3DI2SYiQjuIqHshU7U2LT/i669N0/IjChmb4aKLUjBQSD7ztxzqW9SoNbeOFg47hL1a32WnIiIMFFxKGZvxSouJSptKM+BAb9Lju0M553BkAKIpxdUFj6GiQztWjC+0ONCb28IZbBylFA0/otoKODNepZixqHshArTDmKYfMl3zONSXXzQVubaBgeBHMQ0vZHKhTTuMaQeJFoEoXNOk7ccMljJkbJPBosO3x6qYhkG1FSxWTNUa9eaghYNGcxtxrCQJrtL0Ga+0eW2qzmDRZajoLmoRfXmHaCDPa1N1pmo+vTmH+UZA0Q3oydmrvMP2IIoV9XZItR1QbQfEMfhRjGsb1L2QvGMxVfeQdNvRcoa+vE2skv4IV6sB7SBejNgyRcg6JuWsQ941yTomOScJwgijmFI28eU8eFB4aaJKrCDrmFqj3kS0cNBo7gDlnEPBtZhYaDNV9VhoBewvZxfzJQaLLkopXr5ao9LyyLv2Dau3rtfBe7vww5hqerfe8EKCKMYLYgwDQHAtg3tHS3zr4jxKKVzLpOmFzNZ93tTTwyuTNcIoib5zLIOsbdKfd8i71qIwyNrmdc2R3nakj8fPTCISkHetNMGwsGfKaN9JdLSSRnOHqbUDrlRaBKGir+AwUsosXgAnFlqcGa8SK0XesRjpyXC0K7y128HbHZxwpy6OLT9KBUJAw4toBRFRHGOIYIpBxjYQSQoPmiJEccyl2SZnJmrMNz3yrsXdQwX29SRVa0tZm5xjkrMTgbCWIpLbRTjuBnS0kkazjShmbE4MWUzW2szUfKqtJKKpJ2sz2pMljGLOjNeYb/qYplDIWIvhrd2lOzr76iy/HRdIpRR1L6SadkOrt0MafohSSWhpzrawzUQQiEA7DPFChRfEtMMoNQ9ZfM/dA5Rz9mKhwZxj4lrGhnI6tP/tzqCFg0azBRiGMNqTCIQr8y0uzTYT4VDOcKA3RxTDi+NVxudbWIZQdG2yjrmYEFn3Qvwwpi81w2xmQmQYxYlAaIXMN33qXkg7iECBaxs4ZhJ9FMWKatsnUgo/VIRRTMY2yVgmfXmHvrxNLg3LXWoe0mx/tHDQaLaQnJOYWabrHlNVj9pkwGhPlsP9OYIo5sUrVS7ONMjaJm8YKS2W7vDDmKYfLb6+1Yb3XhhRa4dUmj6zdZ+mn7yHSFJZ1hAhUDELrQBBCOI4zc0wyToW+8sO5Zy9LvOQZnujhYNGs8WICEPFzGLY65X5FpWmz4HeLH4U88KVBc5M1Ci69mJCZBAlF+5bCd9s+ol2MF3zmGt6tPwIL4wQhFgpRCCMQCRZ5toGOcdK+ojk7NRfYKV+Bq0V7Da0cNBotglJYb4Ccw2fiYUW56YbDBZcjg8V+PbYAs9enud7Twzx3nuH+duXJ5mqeuwrZ9ccvhnHirofMlf3may2qbYCKm2fIEwqmgqJucuQJNIoa1vk8xZ9eXvRfJVzLG0e2iNo4aDRbDP68g7FjMVEpc1UzSPnWBztz/HSZJ2nLszyrnuGeNeJIZp+xD0jxZvuK4wSU9DVapvJapuZukelGRBECkHhmgaubVF0E2dxOWfTn3fozTtkHRPX0tnGexUtHDSabYhtGhzqz7HQChivtCjlHAbyDi9PVilmbfaXs9zIktMOkpLWl2ebvD7XYKbu0woilFLkbItS1mKo4DBQcujLOfTlHXqyjjYPaa5DCweNZhvTk7UpuBZXq23iWDFVa/P3T/833lL/D+zjZegD9cBHqAz/OOemGpybqTGx0E7qFUWKnJNoA4f784z2uAwUXPryDoWMrc1DmpuihYNGs80xDWF/OUs5a2NNfYkvPPc3/HX7Lo67Js822ly58DjzfS5x8QQ526KcTxzXh/pyDPdk6M3ZZOzt91PXyWzbm+33jdFoNCuSdy3edOV/wy1FfKL5w/z5wiP0Ww1KZoNB9fccOHqSkXyWYs7GMQ2afsTYfIurC20sU7AMA9sUIElYM0QQktLYki5L/pLlRprYJgDXbX9tPZ0x123ftS9k+XqRNZfy1mwdWjhoNDsIo/U6b8gpfnboS5xtH+C4exnLhFBZBId+izCK8aOYMFYEkaLlRyhAkfxTKExDMA0wJCl1YRmCiGClZqbOa9MApW6P6ekb52bww5goVmll1tub6a1ZP6sKBxHJAF8F3HT7P1NK/Xa67leAXwZC4K+UUv9yydh7gD/tWnQM+FdKqf9bRP534B8D0+m6/00p9dfpuA8DP0/SJvRXlVJf2vAMNZrdRO4Q7foEpgHvKT9Lr1VLlx+GI32LmymlCGNFGCnCOE4fr38exTFBlISxrlRirVMiwzI6tZIMTINUuBgYRlI91TDAEgMxkvwIlUojhSJWybGodH+d55FSDBZdgEXfh259u71Yi+bgAY8opeoiYgNfF5EvAFng/cADSilPRIaWDlRKvQI8CJD2oL4CfLZrk/9LKfX73WNE5F6S3tL3AfuAvxGRE0qpaN2z02h2Gyc/QuPv/lcAckYrWWbm4ORHrttMRLBNIWkdsXo4ahQrgii5k18UKHFHwCTP/SgmDGLieOV9iCRRVolAkUVTVvIoWKaRPBrCkf4c7SC+rvXtZmR6azaPVYWDSsq21tOXdvqngF8EPqqU8tLtplbZ1buBc0qp11fZ7v3An6T7vSAirwFvB76x2rFqNLueoz9Do2pgv/j/xeV8ojGc/Agc/Zlb2m2iDawtpyGOO5qHIohjoih9XBQsiaBpBTfWSjK2yXfGkkZAhYyFF0S0gpjvu2eQuYZ/TaAYiUDZDV3xdhpr8jmkd/3PAHcD/1Yp9aSInAAeFpGPAG3gXyilvnWT3XwA+PSSZb8sIv8T8DTwvyql5oH9wDe7thlLly09pg8BHwI4dOjQWqah0exoOtE9L06+heFjf0zPG7bGeWsYgpNerLNr0Eo6mkcYq0VBMlRy6cnavDheZbbhU8xYvGl/DwrhynxrhfcEy0i0EtuU9DF9bRiYi8Ik0VBWQ0dKrc6ahENq0nlQRMrAZ0Xk/nRsL/AQ8DbgMRE5plZoECEiDvAjwIe7Fv874F+TaCH/GvgD4OdIgyOWHsIKx/Rx4OOQ9HNYyzw0mp1KJ7onaxuUMhYqZsdE91imwUqJ1sOlDG853Hfdso5WsmjW6vKZdExffnjNV7ISInSZtoxFE1dHkMzUPL7y6jQ9WR0pdTPWFa2klKqIyBPAoyR39J9JhcFTIhIDA1xzMHfzg8CzSqnJrn0tPheRfw98Pn05BhzsGnsAGF/PcWo0u41OHwfHMmiHMQMFl3YQ7brono5W4rD63b9SiYBY6htZ6oRPhEm8aN76+3MzeEEEAqWsoyOlbsBaopUGgSAVDFngPcDvkfghHgGeSE1MDjBzg938NEtMSiIyqpSaSF/+KPBC+vwvgT8Wkf+TxCF9HHhqXbPSaHYZnT4OhiR9ICC5O97L0T0iqdN7jeWf4tRH8szrcwyX3MVoLNCRUiuxFs1hFPhU6ncwgMeUUp9PTUWfFJEXAB/4oFJKicg+4BNKqfcBiEgOeC/wT5bs92Mi8iCJyehiZ71S6kUReQw4QxIi+0s6Ukmz1+n0bdDRPRvHMATXMBnpydDyI/1ZroLuIa3R7AC2unf0bkJ/lte4WQ9p3a5Jo9kBjJazvPfeYbKOyUzdI+uYe/Jithnoz3Jt6PIZGs0OYbSc1RewTUJ/lqujNQeNRqPRLEMLB41Go9EsQwsHjUaj0SxDCweNRqPRLEMLB41Go9EsQwsHjUaj0SxDCweNRqPRLEMLB41Go9EsQwsHjUaj0SxDCweNRqPRLEMLB41Go9EsQwsHjUaj0SxDCweNRqPRLEMLB41Go9EsY1XhICIZEXlKRE6LyIsi8jtd635FRF5Jl39shbH3iMjzXX9VEfln6bp/IyIvi8i3ReSzIlJOlx8RkVbXmD/avOlqNBqNZi2spZ+DBzyilKqLiA18XUS+AGSB9wMPKKU8ERlaOlAp9QrwIEDaZvQK8Nl09ePAh5VSoYj8HvBh4NfTdeeUUg9ufFoajUajuRVW1RxUQj19aad/CvhF4KNKKS/dbmqVXb2b5KL/err9l5VSYbrum8CBDRy/RqPRaG4Da/I5iIgpIs8DU8DjSqkngRPAwyLypIh8RUTetspuPgB8+gbrfg74QtfroyLyXLrfh9dyjBqNRqPZPNbUJlQpFQEPpn6Bz4rI/enYXuAh4G3AYyJyTCmllo4XEQf4ERLT0dJ1vwmEwH9OF00Ah5RSsyLyVuBzInKfUqq6ZNyHgA8BHDp0aC3T0Gg0Gs0aWVe0klKqAjwBPAqMAZ9JzU5PATEwcIOhPwg8q5Sa7F4oIh8Efgj4mY5QUUp5SqnZ9PkzwDkSLWXpsXxcKXVKKXVqcHBwPdPQaDQazSqsJVppsCuSKAu8B3gZ+BzwSLr8BOAAMzfYzU+zxKQkIo+SOKB/RCnVXPJ+Zvr8GHAcOL+eSWk0Go3m1liLWWkU+FR6wTaAx5RSn09NRZ8UkRcAH/igUkqJyD7gE0qp9wGISA54L/BPluz3DwEXeFxEAL6plPoF4J3A74pICETALyil5m55phqNRqNZM7KCi2DHcerUKfX0009v9WFoNBrNjkJEnlFKnVpx3W4QDiIyDby+1cdB4nO5kWltN7Cb57eb5wZ6fjud2zW/w0qpFZ22u0I4bBdE5OkbSeHdwG6e326eG+j57XS2Yn66tpJGo9FolqGFg0aj0WiWoYXD5vLxrT6A28xunt9unhvo+e107vj8tM9Bo9FoNMvQmoNGo9FolqGFg0aj0WiWoYXDKojIn3Y1HrqYVqftrHtARL6RNjv6johkVhi/YlOjdN2HReS1tGHSD9yZGS07vlud30+k62MROdW1fFs0bbpd80vX7Ybz1ycij4vI2fSxN12+W87fivNL123p+bvR3Nb62YvIyXT+3xGR/yoipfWMXxWllP5b4x/wB8C/Sp9bwLeBk+nrfsBcYcz3A1b6/PeA30uf3wucJikhcpSkwOCy8Ttgfm8E7iEpyHiqa/kR4IWtPme3cX675fx9DPiN9PlvdH0/d8v5u9H8ttX5WzK3NX32wLeAd6XPfw7415t57rTmsEYkKQD1k1wrIPj9wLeVUqcBlFKzKiltfh3qxk2N3g/8iUqq0F4AXgPefjvncDNuYX4vqaTj37bmNsxvV5w/knl8Kn3+KeAf3uZD3RC3YX7b5vytMLe1cg/w1fT548A/2szj0sJh7TwMTCqlzqavTwBKRL4kIs+KyL9cwz66mxrtBy53rRtLl20VmzG/pRyV7dO0abPnt1vO37BSagIgfexu97sbzt+N5redzt/SucHaPvsXSPrkAPwEcHCd42/Kmpr97HZE5G+AkRVW/aZS6i/S50vLjlvAPyBpdNQE/laSIlZ/e4P3WNrUSFbY7LbEFd+J+a3Ampo2bQZbND99/jaJ3Xz+Nji3tX72Pwf8f0TkXwF/SVIdez3jb4oWDoBS6j03Wy8iFvBjwFu7Fo8BX1FKzaTb/DXwFmDZl1OuNTV6t0qNgun4bkl/ABjf6Bxuxu2e3w3e0wM6/cWfEZFO06ZNL5+7FfNj95y/SREZVUpNiMgoSSvg3XT+Vpwfd+j8bWRua/3slVIvk5jXOj11/of1jF8NbVZaG+8BXlZKjXUt+xLwgIjk0hP8LuDM0oFyg6ZGJJL+AyLiishRkqZGT922GdycDc/vRsj2atq06fNj95y/vwQ+mD7/IPAXsKvO34rzY/ucv2VzW+tnLyJD6aMB/BbwR+sZvyq36tHeC3/AfyRpOrR0+c8CL5LY/j7WtfwTpJEtJI6uy8Dz6d8fdW33myRREq8AP7hD5/ejJHdhHjAJfCld/o/SsaeBZ4Ef3k3z20Xnr5/kbvts+ti3y87fivPbLudvpbnd7LNfMrd/Crya/n2UaxUvNuXc6fIZGo1Go1mGNitpNBqNZhlaOGg0Go1mGVo4aDQajWYZWjhoNBqNZhlaOGg0Go1mGVo4aDQajWYZWjhoNBqNZhn/f5JpevDpfIv1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "coords = []\n",
    "for r in range(assignment.shape[0]):\n",
    "    coords.append([(assignment.Base_Longitude[r],assignment.Base_Latitude[r]),(assignment.Inc_Longitude[r],assignment.Inc_Latitude[r])])\n",
    "# plot the line segments, indicent points, and base station points of the final network\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlim(-76.23, -75.93)\n",
    "ax.set_ylim(36.72, 36.93)\n",
    "lc = mc.LineCollection(coords, alpha = .2)\n",
    "ax.add_collection(lc)\n",
    "plt.scatter(data = ohca_df, x = 'Longitude', y = 'Latitude', alpha = .3)\n",
    "plt.scatter(data = stations_df[stations_df['Base'].isin(bases_selected.index.tolist())], x = 'Longitude', y = 'Latitude', alpha = 1, c = 'orange') \n",
    "plt.scatter(data = stations_df[~stations_df['Base'].isin(bases_selected.index.tolist())], x = 'Longitude', y = 'Latitude', alpha = 0.2, c = 'orange') \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9585e22c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Freeing default Gurobi environment\n"
     ]
    }
   ],
   "source": [
    "m.dispose()\n",
    "gp.disposeDefaultEnv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4e2add8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
